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NOTATION 


A list of the symbols used throughout this document and their definitions is 
provided below for convenience. 


Roman Symbols 


a. .. speed of sound 
cp . . . specific heat at constant pressure 
Cy ... specific heat at constant volume 
e . . . internal energy 

i . .. z index of numerical solution 

j . .. r index of numerical solution 

k . .. 0 index of numerical solution 

n. . . time step index of numerical solution or rotational speed (rpm) 
n... outward unit normal vector 
p . . . pressure 

r ... radius or radial coordinate 
t ... time 
v ... velocity 
z ... axial coordinate 
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A... surface area 

9 

Cp ... pressure coefficient (Cp = p — poo/2poo«oo) 

Cp ... power coefficient ( Cp = P/pn^D^) 

Cf . . . thrust coefficient (C^ = T/pn^D^) 

CFL . . . Courant-Freidrichs-Levy 

D ... dissipation flux vector or diameter 

DELT . . . adjacent cell grid spacing 

DELTLE . . . adjacent cell grid spacing at leading edge 

DELTTE . . . adjacent cell grid spacing at trailing edge 

F . . . flux vector in z direction 

G . .. flux vector in r direction 

H . .. flux vector in 6 direction 
Hi . . . total enthalpy 

J . . . advance ratio (J = U/nD) 

K . . . source term flux vector 
L .. . length 

N . .. number of propeller blades 
0 ... orthogonality 

P ... point distribution function or power 
Q ... vector of dependent variables 

R. .. gas constant or residual or maximum radius 
RAT . . . maximum ratio of adjacent cell grid spacings 

S . .. arc length or pertaining to surface area normal 
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T . . . temperature or torque 

U . .. flight velocity 

V . .. volume 


Greek Symbols 

a . . . time-stepping factor 
fi . . . 3/4 radius propfan blade setting angle 
t . . . modified second-order damping coefficient 
. . . modified fourth-order damping coefficient 
p ... density 

k l . . . second order damping coefficient 
. . . fourth order damping coefficient 

7 . . . specific heat ratio 

5 .. . spatial second-order centred difference operator 
A . . . blockage factor 

T] ... radial transformed variable for C-grid generation 
£ . . . axial transformed variable for C-grid generation 
v ... damping factor 

T . . . relative pressure gradient coefficient 
A... increment of change 

Special Symbols 

V . . . spatial vector gradient operator 
A... spatial forward difference operator 
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V . . . spatial backward difference operator 

Superscripts 

[ ] . . . averaged variable 

[~] . . . dimensional variable 

[ ] . . . implicitly smoothed variable 

[ J . . . vector variable 

[ ]* . . . intermediate variable 

[ ] n . . . time step index of variable 

Subscripts 

[ ]j j . . . grid point index of variable 
[ )max • • • maximum value 
[ Imin • • • minimum value 
[ )p . . . related to pressure 
[ ]p 5 . . . pressure (high pressure) surface 
[ ] 55 . . . suction (low pressure) surface 
[ . . . total quantity 

[ ] z . . . derivative or value with respect to z 
[ ]r . . . derivative or value with respect to r 
[ . . . derivative or value with respect to 6 

[ ]qo . . . freestream value 
[ ] re y . . . reference value 



1. SUMMARY 


The time-dependent three-dimensional Euler equations of gas dynamics have 
beem solved numerically to study the steady compressible transonic flow about ducted 
propfan propulsion systems. Aerodynamic calculations were based on a four-stage 
Runge-Kutta time-marching finite volume solution technique with added numerical 
dissipation. An implicit residual smoothing operator was used to aid convergence. 
Two calculation grids were employed in this study. The first grid utilized an H-type 
mesh network with a branch cut opening to represent the axisymmetric cowl. The 
second grid utilized a multiple-block mesh system with a C-type grid about the cowl. 
The individual blocks were numerically coupled in the Euler solver. Grid systems were 
generated by a combined algebraic/elliptic algorithm developed specifically for ducted 
propfans. Numerical calculations were initially performed for unducted propfans to 
verify the accuracy of the three-dimensional Euler formulation. The Euler analyses 
were then applied for the calculation of ducted propfan flows, and predicted results 
were compared with experimental data for two cases. The three-dimensional Euler 
analyses displayed exceptional accuracy, although certain parameters were observed 
to be very sensitive to geometric deflections. Both solution schemes were found to 
be very robust and demonstrated nearly equal efficiency and accuracy, although it 
was observed that the multiple-block C-grid formulation provided somewhat better 
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resolution of the cowl leading edge region. 
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2. INTRODUCTION 


Modern turbopropeller (propfan) aircraft propulsion systems have demonstrated 
the potential for significant fuel savings over existing high bypass turbofan engines. 
By incorporating innovative airfoil design concepts utilizing a high degree of blade 
lean, sweep, and twist, propfans can perform effectively at relative Mach numbers 
extending well into the transonic regime. A number of aerodynamic phenomena com- 
mon in propfan installations are illustrated on Fig. 2.1. The process of designing 
advanced propfans requires a complicated iterative interaction between the disci- 
plines of aerodynamics, structural analysis, and aeroacoustics. Reliable attainment 
of performance goals necessitates the use of proven aerodynamic analyses for the 
highly three-dimensional blades. In this respect, the realization of the performance 
improvements available in advanced propfan designs is directly related to advances 
in computational fluid dynamics. 

Several researchers have recently published advanced numerical analyses for prop- 
fan flows. Barton et al. [1] and Usab et al. [2] presented results from several three- 
dimensional Euler codes for single-rotation propellers which were verified through 
comparison with experimental data. Celestina et al. [3] performed a time-averaged 
solution of three-dimensional counter-rotating propeller flowflelds using an average- 
passage system of equations. Whitfield et al. [4] performed calculations for unsteady 
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three-dimensional single-rotation and counter-rotation propfans using a multiple- 
block grid formulation. Saito et al. [5], Matsuo et al. [6], and Kobayawa and Hatano 
[7] have all demonstrated the application of the three-dimensional Navier-Stokes equa- 
tions for the analysis of high speed propfan designs. Many of the aerodynamic features 
relevant to propfan flowfields have been successfully captured by these and other anal- 
yses, including blade leading edge and tip vortices, transonic/shock flow phenomena, 
and counter-rotating blade row unsteady interaction. 

In light of the desire to optimize propfan designs, the concept of a ducted prop- 
fan configuration (very high bypass fan, see Fig. 2.2) has recently been explored. By 
ducting the propulsor, a significant increase in thrust and efficiency can be obtained, 
particularly at static, or low forward flight velocities. This thrust improvement could 
significantly impact future designs by permitting reduced diameters and underwing 
installation. Ducted propellers may also have advantages related to community noise 
levels and enhanced public acceptance. The concept of a ducted propeller has existed 
for more than fifty years, although research was initially directed at low-speed appli- 
cations due to the drag penalty associated with the cowl at higher flight velocities. 
Early tests of ducted propellers are given in a review by Sacks and Burnell [8]. A 
more recent discussion on the technical issues relating to ducted propfans was given 
by Groeneweg and Bober [9]. In order to assess the relative aerodynamic merits 
of advanced ducted propfan design configurations, detailed computational analyses 
similar to those described above for unducted propulsors are required. 

This document contains the Final Report and User’s Manual for the 3D Eu- 
ler Ducted Propfan Analysis codes developed by the Allison Gas Turbine Division 



PROPFAN AERODYNAMIC CHARACTERISTICS 


5 



6 


ORIGINAL PAGE 
COLOR PHOTOGRAPH 



Figure 2.2: Ducted propfan (very high bypass fan) geometry 
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of the General Motors Corporation under Task 1 of NASA Contract NAS3-25270. 
This study was directed at the development of a three-dimensional Euler analysis for 
high-speed ducted propfans. A finite volume multiple-block four-stage Runge-Kutta 
numerical algorithm was utilized to predict the aerodynamics of ducted fan flows. Of 
particular interest was the ability to accurately predict the aerodynamics of the cowl. 
Two different grid arrangements were tested. The first system involved a sheared 
H-type grid which utilized a branch cut opening for the cowl. The second system was 
based on a numerically coupled multiple-block grid arrangement with a body-centered 
C-type grid about the cowl. Predicted results using both grid systems were compared 
with experimental data for two cases: a single-rotation low speed ducted propeller, 
and a high-speed 1.15 pressure ratio fan (both single-rotation and counter-rotation 
arrangements). This program is considered the initial step in developing a complete 
aerodynamic model of the flowfield through a ducted fan and the external flow over 
the engine nacelle and fan cowl. 

To predict the flow about a ducted propfan using the analyses described in this 
document, the necessary sequence is: 

1. Generate a numerical grid for the domain of interest. 

2. Run the Euler code to predict the steady aerodynamics. 

3. Process the results as needed. 

Separate sections are provided in the chapters which follow to describe the basis 
and operation of the codes used in the first two steps. A theoretical development of 
the grid generation algorithms for both the H-grid and C-grid architectures is given in 
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Chapter 2. The 3D time-marching Euler analysis is detailed in Chapter 3. In Chapter 
4, a summary of the predicted results and verification studies performed to validate 
the accuracy of the analysis is presented. A summary of the conclusion of this study 
is given in Chapter 5. A user’s manual for the grid generation code is provided in 
Appendix A. The operation of both the H-grid and multiple-block C-grid 3D Euler 
computer codes is outlined in Appendix B. 

It is worthwhile mentioning that most of the plotting and graphical postprocess- 
ing of the solutions are performed on graphics workstations. Presently, the PL0T3D 
[10] graphics software developed at NASA Ames Research Center is being extensively 
used for this purpose, and plot output has been tailored for this software. 
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3. GRID GENERATION ALGORITHM 


In this chapter, the numerical algorithm forming the basis of the ducted propfan 
analysis grid generation schemes is described. The geometry and computational do- 
main are briefly described in the first section below. During the course of this study, 
two types of grids were tested for the numerical calculations. The grid systems are 
divided into two categories: H-grid, and multiple-block C-grid. Separate sections are 
provided to describe each category. 

3.1 Computational Domain 

The problem of interest is a ducted propfan geometry in steady operation as 
shown in Fig 2.2. Geometric parameters are expressed in a cylindrical coordinate 
system. The axial coordinate lies along the propeller axis of rotation, the radial 
coordinate is perpendicular to the axis of rotation, and the circumferential coordinate 
sweeps in the counter-clockwise direction when viewed down the axis of rotation. This 
coordinate system is illustrated in Fig. 3.1. The hub and duct contours are assumed 
to be axisymmetric surfaces (no circumferential variation). It is further assumed that 
the airfoil tips are fully enclosed within the duct (although this limitation could be 
easily relaxed). For any number of blades, under the assumptions of steady flow and 
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a periodic geometry, the computational domain may be reduced to a single blade 
passage, circumferentially, as shown in Fig. 3.1. For counter-rotating configurations, 
a single blade passage is used for each blade row, and the multiple passages are then 
numerically coupled in the flow solver. The problem is further simplified by fixing the 
computational domain to the rotating blade elements. The resulting flow predictions 
are therefore based on the steady flow relative to the rotating blade. 

In order to implement a finite volume numerical solution of the steady aero- 
dynamics for this computational domain, the region of interest must be subdivided 
into a finite number of smaller elements within the overall region of interest. To 
enhance the quality of the numerical solution, we further specify that these subdivi- 
sions be constructed in a relatively smooth and orderly manner with some arbitrary 
limitation on both the relative spacing between cells and the total number of cells 
in the computational domain. Each of these subdivisions will be used to define an 
elemental computational cell which will form the geometric basis of the finite volume 
aerodynamic solver to be described in a later chapter. 

In this study, two methods for constructing the computational cells were devel- 
oped. The first method, referred to as the H-grid, is based on a single-block array of 
computational cells ordered in an H-type fashion. The second method, referred to as 
the multiple-block C-grid, utilizes a coupled system of five blocks each containing an 
array of computational cells. The five-block system is constructed by first subdividing 
the overall computational domain into five subregions, then further subdividing the 
individual subregions. The multiple-block arrangement permits the use of a special 
C-type grid wrapped about the leading edge of the duct to enhance the accuracy of 
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Figure 3.1: Ducted propfan analysis computational domain 
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the aerodynamic calculations in this critical region. 

Both the H-grid and C-grid algorithms are described in the sections below. 

3.2 H-Grid Algorithm 

The H-grid generation scheme is based on a procedure originally developed by 
Mulac [11] for multiple blade row flow calculations using the average-passage [12] 
system of equations. In solving the average-passage system of equations, it is a re- 
quirement that the individual blade row grids have identical axisymmetric projections. 
This same technique was also applied in this study. The grid generation technique 
applies equally to single- rotation or multiple- rotation geometries. The approach here 
is to first construct a two-dimensional axisymmetric mesh in the z,r plane, followed 
by separate constructions for each of the individual blade rows in the circumferential 
direction to determine the final three-dimensional meshes. 

The two-dimensional meshes are constructed by dividing the axisymmetric plane 
into several subregions as shown in Fig. 3.2 for an unducted propfan and Fig. 3.3 for a 
ducted propfan. The regions of interest include the inlet, mid-inlet, blade, mid-exit, 
exit, between blade, tip gap, and outer flow regions. The mid-inlet and mid-exit 
regions result from the assumption that the blade leading and trailing edges are fully 
contained within the axial span of the duct. These regions are individually gridded 
in the two-dimensional plane to satisfy the conditions of common grid points along 
region interface boundaries, and a suitable distribution of points within each region. 
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Unducted Propfan H-Grid Generation Subregions 



Between Blade Region 


Figure 3.2: Axisymmetric plane projection grid generation subregions for an un- 
ducted propeller 
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Ducted Propfan H-Grid Generation Subregions 



Between Blade Region 


Figure 3.3: Axisymmetric plane projection grid generation subregions for a ducted 

propeller 
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The various regions differed in their individual construction. The distribution of 
points in any region is given axially and radially by one of three point distribution 
routines. These routines are described in the paragraphs below. All interpolations of 
coordinates were performed using a cubic spline interpolation routine. 

The distribution of the points on the axisymmetric blade plane in both the radial 
and axial direction is determined by a simple geometric progression radiating from a 
symmetric centerline referred to as packing algorithm #1 described below. 


Packing algorithm =& 1 


Given: 

D total projected length of cubic spline curve 
RAT maximum ratio of adjacent cell projected lengths 
M number of points distributed across D (must be odd) 
the initial cell length of the progression is calculated as 

The point distribution is then given by 


(3.1) 


P( i) = E (DELT)(RAT)3 ; i = 1, (M/2) + 1 (3.2) 

j = 1 


Packing algorithm #2 


The second construction utilizes a slightly different packing algorithm. In this 
case, D, RAT , and DELT are specified along with N, the total number of parallel 
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spline curves, leaving M to be determined. The value of M is calculated iteratively 
until the the following conditions are satisfied: 

D \ 


(,DELT ) ma x ^ RAT ? 

yy m .„=2(^r 


where 


f— ) 

\deltj max 

D \ 


max 


[— 

\DELT 


);=!,« 


( DELT ) min m ‘” {dELt) j=l,l 


( 3 . 3 ) 


MELT} min \DELT J 7 = 1 , N 

These conditions guarantee that neither the adjacent cell ratios, nor their in- 
verses, exceed RAT. The resulting point distribution is then given by 


P(i) = 53 (DELT)(RATy ,* = 1 ,M 
j = 1 


( 3 . 4 ) 


This approach is used axially in the inlet and exit regions, and radially in the 
outer boundary region. 


Packing algorithm #3 


The third packing algorithm is used to determine the axial point distribution 
between blade rows. 

Given: 

D total projected length of cubic spline curve 
RAT maximum ratio of adjacent cell projected lengths 
DELTTE initial cell width (forward blade row trailing edge) 



17 


DELTLE initial cell width (aft blade row leading edge) 

N total number of parallel cubic spline curves 
the following values are then calculated: 

dte -{d elt D tb7dbltle )^ < 3 - 5 > 

dlb =( deltteTd E eltle )( d > < 36 > 

The value of M is determined iteratively until the following conditions are satisfied: 

/ n \ M 

(mif (3.n 
(: oEuf < 3 - 8 > 

where: 

( D \ ( DTE \ ( DLE \ . mr , v 

\DELTJmax ~ m<1X \DELTTEJ \DELTLEJ ' 3 ~ N . 

/ D \ . ( DTE \ ( DLE \ „ 

KDELTJmin ~ \DELTTe) \DELTLEJ ~ lfN ( 31 °) 

Again, this ensures that neither the adjacent cell ratios, nor their inverses exceed 

RAT. The point distribution is then given by: 

t 

P(i) = £ ( DELTTE){RATy , » = 1, M (3.11) 

j=l 

P(i) = £ ( DELTLE)(RAT)J , * = 1, M (3.12) 

3 - 1 

The blade region is generated initially. This region is defined by the axisymmetric 
projection of the blade hub and tip boundaries, and the blade leading and trailing 
edges. The number of points defining the blade region axially and radially is specified. 
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The point distributions are determined by packing algorithm #1. Once the blade 
points are determined, the full inlet region pictured in Fig. 3.2 is computed next using 
packing algorithm #2, regardless of whether a ducted or unducted configuration is 
being generated. Once this region has been discretized, a simple linear interpolation is 
performed to determine the coordinates of the radial grid line which exactly intersects 
the duct leading edge. This grid line then determines the boundary between the mid- 
inlet region and the inlet region in Fig. 3.3. Once this boundary is determined, the 
upper boundary of the mid-inlet region is constructed by extrapolating the blade tip 
grid line axially, and maintaining the known tip gap spacing. Packing algorithm #1 
is then used to distribute points axially in exactly the same manner in the mid-inlet 
region as the blade region, with the additional constraint that the grid lines must 
match at the region interface. Once this region is constructed, the inlet region can be 
constructed using packing algorithm #2. The point distributions for additional blade 
rows (if necessary) are then determined using packing algorithm #1. For multiple 
blade rows, the next step is to calculate the region between blades using packing 
algorithm #3. The mid-exit and exit regions are constructed in the same manner as 
the mid-inlet and inlet regions by locating the grid line which exactly intersects the 
cowl trailing edge. The tip gap region is constructed algebraically with a specified 
number of equally spaced points in the radial direction between the blade tip grid 
line and the cowl lower surface. Upstream and downstream of the cowl, the cowl 
surface grid line is extended axially from the leading and trailing edges. It should be 
noted that this may impose some limitation on the maximum diameter of the hub 
upstream and downstream of the cowl. Finally, the radial distribution of points in 
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the outer region is determined from packing algorithm #2. The grid lines which form 
the upper and lower surface of the cowl overlap upstream and downstream of the 
cowl. A duplicate grid line is placed on top of the cowl lower surface grid line to 
provide a natural storage location for the phantom points used in the application of 
the cowl surface boundary conditions (see Section 4.6). Once the axisymmetric plane 
gnd has been calculated, the full three-dimensional grid is constructed individually, 
for each blade row. The blade twist and thickness distributions are used to determine 
the circumferential coordinates for the blade regions. The remaining circumferential 
coordinates are based on utilizing packing algorithm #1 in the blade-to-blade plane. 
The circumferential variation of grid points upstream and downstream of the blade 
regions is adjusted to provide a smooth transition such that the grid lines become 
parallel to the axis of rotation away from the blade region. 

A sample E-grid for a ducted propfan based on this grid generation technique 
is given in Fig. 3.4. The tip gap has been enlarged for this geometry in order to 
illustrate the tip gap grid. The emphasis in this development, as in Mulac [11], was 
to maintain a reasonable grid quality, rather than a specific number of grid points. 
Due to the algebraic nature of the scheme, the construction is performed relatively 
easily and inexpensively. Unfortunately, this approach can have drawbacks under 
certain conditions. The cubic spline interpolation routine is subject to errors when 
a boundary or grid line has a discontinuous or sudden change in slope (particularly 
troublesome spots are the spinner or cowl leading edges). In addition, blades or cowls 
with a large thickness or blunt leading edge may not be well represented by this 
scheme when a relatively low number of grid points is used. Finally, the construction 
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may occasionally have limits based on the number of grid points in a particular region 
due to the inability to satisfy one or more of the adjacent cell ratio constraints in the 
packing algorithms. 


3.3 C-Grid Algorithm 

The C-grid generation scheme utilizes a number of concepts developed for the 
H-grid scheme. Again, the mesh is generated by first constructing a two-dimensional 
grid in the axisymmetric plane. In this case, a slightly different subdivision approach 
is utilized. The individual subregions and the five-block structure for the C-grid 
approach are listed and displayed in Fig. 3.5. Several of the subdomains are similar 
to those used in the H-grid construction, so the extension of the algorithm to the 
C-grid is obvious. In the case of the H-grid, the individual subregions are eventually 
combined into a single grid for use with the 3D H-grid Euler analysis. In the multiple- 
block C-grid generation, the subdomains are left divided and the individual flow 
solutions are numerically coupled. The blade, mid-inlet, inlet, mid-exit, exit, and 
outer regions are all constructed in basically the same manner for the C-grid as they 
were for the H-grid. In this case, the axial locations of the upstream boundary of 
the mid-inlet and the downstream boundary of the mid-exit regions are specified as 
a fraction of the cowl axial chord upstream and downstream of the cowl leading and 
trailing edges, respectively. The tip gap is extended upstream and downstream of the 
cowl, so there are no overlapping grid lines (resulting in the need for blocks 2 and 4), 
and the outer boundary extends radially Horn a grid line which is displaced from the 
cowl outer surface by a distance equal to the inner surface tip gap spacing. The near- 
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Figure 3.4: Sample H-grid for a ducted propfan 
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cowl region is then discretized by the C-gnd as shown in Fig. 3.5. The construction of 
blocks 1,2,4, and 5 follows directly from the development of the H-grid, so no further 
discussion is presented here. Instead, the underlying algorithm for the C-grid (block 
3) is presented in the paragraphs below. 

Grid points for the C-grid block surrounding the cowl are determined in a two- 
step procedure. In the first step, the inner and outer boundary points are specified. 
In the second step, the interior points are determined. Each step has several options 
each of which is described below. 

The exterior boundary points of the C-grid are determined by the surrounding 
four blocks. This is based on the requirement of coincident points along all block 
boundaries. The number of points along the inlet plane is directly specified, and 
must be an odd number. The number of points from the cowl surface to the outer 
boundary of the C-grid must also be specified, which fixes the number of points along 
the C grid block exit plane. The cowl surface points axe determined through one 
of two procedures. The first method performs a simple interpolation of the input 
coordinates based on arc length around the airfoil such that the relative arc length 
distributions of the airfoil surface points and the corresponding outer boundary points 
are the same. Unfortunately, this can often lead to highly skewed mesh lines near the 
airfoil surface. To circumvent this problem, an option was provided which allows the 
cowl surface grid points to “float” along the contour of the airfoil in a manner that 
would impose orthogonality at the cowl surface from a line extending from the cowl 
point to the corresponding outer boundary point. The movement of the airfoil grid 
points is controlled by a secant iteration procedure which optimizes the orthogonality 
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Ducted Propfan Multiple-Block C-Grid Subregions 



Blade Region 


Figure 3.5: Axisymmetric plane projection subregions for multiple-block C-grid gen- 

eration 
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of the surface grid point location as a function of arc length along the airfoil while 
maintaining a smooth transition through neighboring surface grid points. The new 
airfoil coordinates are determined from the updated value of arc length through a 
linear interpolation of the arc length and airfoil coordinates originally specified. The 
resulting grid thus possesses more desirable orthogonality characteristics along the 
airfoil contour, which is desirable in terms of solution accuracy. 

The secant iteration procedure described above for the cowl surface grid points 
is expressed as: 


s k + 1 _ s k , (0.0 -oht^-sf- 1 ) 

* i+ (oh or 1 ) 


( 3 . 13 ) 


where is arc length measured clockwise around the cowl from a fixed reference 
location (1) to the point (t), k is the secant iteration count, rind O t - is the measure of 
nonorthogonality. The orthogonality measure was based on the change of arc length 
Sgrid between the outer point and the corresponding inner point. The optimized 
point occurs when Oj = dSg r ^/dS^ = 0. 

The starting values for the secant iteration are determined by using an initial 
point distribution obtained from the relative arc length cowl surface grid point inter- 
polation routine (method one, above). 

In order to avoid overlapping grid lines and to maintain stability, the new surface 
grid point locations were never allowed to migrate more than one third of the distance 
between the previous surface location and that of the neighboring grid points. 


Interior grid points are determined through one of two methods. The first method 
is a simple algebraic interpolation of coordinates between the inner and outer C-grid 
boundary point distributions. In the second method, the interior points are generated 
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through the numerical solution of a set of elliptic equations controlling a weighted 
distribution of grid smoothness, orthogonality, and grid point density based on a 
variational formulation originally developed by Brackbill and Saltzman [13]. A brief 
description of this scheme is given below. 

For a two-dimensional grid, the following integral expressions may be derived to 
evaluate critical aspects of the overall grid quality in physical space: 


Grid smoothness: 



h = f J [(Vf) 2 + (Vr,) 2 \dzdr 

(3.14) 

Grid orthogonality: 

I o = j j [(VO • (Vi))]dzdr 

(3.15) 

Grid point density: 

Iw — j J , r)Jdzdr 

(3.16) 

where: 
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(3.18) 


and where V is the Cartesian gradient vector operator. The term w(z,r) is a user- 
specified function the magnitude of which is proportional to the desired grid point 
density in physical space. 


Obviously, the smoothest possible grid is obtained when I s is minimized, the 
most orthogonal grid is obtained when /<> is minimized, and the grid with the most 
desirable point density is obtained when Iw is minimized. By minimizing a weighted 
sum of these terms, i.e: 


I = Is + Colo + C W I W 


(3.19) 



26 


the constants Co, Cw may be used to control the relative importance of orthogonality 
and point density, respectively, in the overall grid point distribution. By exchang- 
ing dependent and independent variables, and applying the concepts of variational 
calculus for minimizing functions using the Euler-Lagrange equations, the following 
nonlinear coupled set of equations results: 

o 1 dxo 

h z (i + b 2 z ( v + b Z z VV + a l r tt + a 2 r ( v + a Z r VV = - J 2^07 ( 3 - 2 °) 

o 1 dw . 

a\ z X + ^2 Z ( V + a Z z VV + c l r tf + c 2 r (Tj + c 3 r W = ~ J ( 3 - 21 ) 

where the coefficients (aj, 6 j,c t -,i = 1,3) are all functions of the coordinate derivatives 
as: 

a l = a sl+Co< l ol + C‘w a vl a 2 = a s2+^o a o2"t"^^ a t;2 a 3 = a a2+Co a oZ+Cw<*’vZ 

b\ = b sl + C Q b ol -I- Cwb v i b 2 = b s2 + C o b 0 2 + C w b v 2 63 = b s 3 + C 0 b 0 3 -I- Cwb v 3 

C 1 = c sl + CoC 0 1 + Cwc v i C 2 = c s 2 + Coc 0 2 + Cwc v 2 C 3 = c s 3 + C 0 c 0 3 + C w c vZ 

(3.22) 

where: 

a sl = “(««)« a s 2 = 2(oa)/3 a a3 = -(aa )7 

b s i = (bb)a b s2 = -2(bb)f3 b sZ = (66)7 

c 5 l = ( cc ) a c s 2 = - 2 («O 0 c s3 = -(“h 
a ol = Z V T V a o2 = z ( r V + z V r ( a oZ = z £ r ( 
b ol = z v b o2 = 2 ( 2z £ z 7 + 6 o3 = 

c ol = r ? c o2 = 2 ( Z ( Z V + 2r £ r »?) c o3 = r \ 



27 


a vl = ~ z V rr l a v2 ~ z ( r l + z V r { a vZ ~ ~ z ( r i 
b vl = r % b v2 = -2r^r v b vZ = 

C v\ = z »7 ^2 = —2 z £ s Tj C|;3 = 

(aa) = z^ + z v r v ( bb ) = r| + (cc) = z| + 

a = ( x % + yrj)/J^ P — ( x £ x rj + y^yr])/J 2 7 = (x^ + (3.23) 


When Co — 0 and Cto — 0 the more standard Laplace or Poisson grid generation 
schemes result. This system is more complex than the usual Poisson-based grid 
generation schemes [14], but is still solvable using standard relaxation techniques. 
In this case, an iterative successive overrelaxation Gauss-Seidel solution technique is 
applied to solve the finite-difference equations resulting from a second-order central- 
difference approximation of the resulting equations. For example: 


z et « 

« (A£)2 

*i+l,j+l - *i+l,j-l + *.-1,7-1 - *.-1.7+1 
(A£A,) 


(3.24) 

(3.25) 


A sample C-grid mesh network is given in Fig. 3.6 for a ducted propfan geom- 
etry. While this arrangement is not remarkably different from the H-grid network, 
the elimination of the highly skewed grid lines from the cowl leading edge region is 
certainly a desirable quality. 
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Figure 3.6: Sample multiple-block C-grid for a ducted propfan 



4. 3D EULER NUMERICAL ALGORITHM 


In this chapter, the numerical algorithm forming the basis of the ducted prop- 
fan aerodynamic analysis is described. Separate sections are presented for nondi- 
mensionalization, governing equations, boundary conditions, and discretization and 
time-integration schemes. 


4.1 Nondimensionalization 


To simplify the numerical treatment, all variables in the numerical solution are 
nondimensionalized by reference values as follows: 


z — 


'ref 


r = 


P = 


^ref 

V 
Pref 


v z = 


P = 


Vz 

v ref 
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Pref 


Vx — 


v r 


v 9 

Vf) = — — 

v ref v ref 


V ref 


(4.1) 


The reference quantities are defined as follows: 


'ref 


is determined from the maximum diameter of all the blade rows 


P re f is determined from the freestream toted pressure 
P re f is determined from the freestream total density 
a Tt j is determined from the freestream conditions 
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= y/lPref/Pref 

v re j is determined from the freestream acoustic velocity as 



4.2 Governing Equations 


The Euler equations of gas dynamics form the basis for the numerical solution 
procedure. In cylindrical coordinates, the strong conservation form of the inviscid 
equations of gas dynamics may be expressed in a reference frame rotating with the 
propeller as: 


J —(XQ)dV + + \OdAr + X(H — — J XSdV + J XKdV 

(4.2) 


where: 
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outside blade row; 
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(4.6) 


Here, v z ,vr,VQ represent the velocity components in the axial, radial, and circumfer- 
ential directions, respectively. The total energy function, e, is defined as: 


The total enthalpy, H y is related to the total energy by: 


(4.7) 


B = H + ~ (4.8) 

P 

These equations are based on the average-passage system of equations developed by 
Adamczyk [15] for the prediction of time-averaged flows in counter-rotating and mul- 
tiple blade row machines. The system of equations is closed through the specification 
of the blockage (A) and source terms (5). The closure model developed by Adamczyk 
[16] was utilized in this study to define these terms. Further details are given in the 
references cited. 

The governing equations are applied to a generalized finite volume in physical 
space as shown in Fig. 4.1. The cell surface areas dA z , dA r , and <LAq are calculated 
using the cross product of the diagonals of a cell face, and the cell volume is determined 
by a procedure outlined by Hung and Kordulla [17] for generalized nonorthogonal 
cells. 

For convenience, a special operator L{Q ) is defined as: 


L{Q) = — \^FdAz + A GdAr + A (H — rwQ)dAQ^j + j A SdV + J XKdV (4.9) 
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i+1 ,j+1 ,k+1 
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The flux variables F,Q, and H are determined at each grid cell interface by averaging 
the cell-centered dependent variables from the individual finite volumes adjoining the 
interface. The L operator represents the summation of fluxes of mass, momentum, 
and energy for the control volume V under the influence of blockage and source terms. 

4.3 Runge-Kutta Time Integration 

The time-stepping scheme used to advance the discretized equations is a four- 
stage Runge-Kutta integration. The solution proceeds as: 

Ql=Q n - *iA WQ n ) + D(Q n )) 

Q 2 =Q n -a 2 At[L(Q 1 ) + D(Q n )l 
Q i = Q n - a 3 At[L(Q 2 ) + D(Q n )) 

Q 4 = Q n -a i At{L(Q i ) + D(Q n )) 

<? n+1 = Qi (4.10) 

where: 

111 

a l = g q 2 = 4 a 3 = 2 Q 4 = 1 (4.11) 

An artificial dissipation operator (D) has been added to control oscillations in the 
numerical solution. A further discussion on this term is given in the following section. 
Jameson [18] determined that this scheme is stable for all time steps satisfying the 
CFL - related time step limitation 


CFL < 2V2 


(4.12) 
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An acceleration technique known as local time stepping is used to enhance convergence 
to the steady-state solution. Local time stepping utilizes the maximum allowable time 
increment at each point during the course of the solution. While this destroys the 
physical nature of the transient solution, the steady-state solution is unaffected and 
can be obtained more efficiently. 


4.4 Artificial Dissipation 


The Runge-Kutta algorithm described above utilizes an artificial numerical dis- 
sipation term D(Q) to suppress the odd-even point decoupling often observed for 
central-difference formulations. This problem is especially prevalent near shock waves, 
and it has been observed that the formulation of the dissipative term can have a sig- 
nificant influence on the final numerical solution. Jameson [18] demonstrated that a 
dissipative system combining second and fourth difference smoothing terms can effec- 
tively eliminate undesirable numerical oscillations without destroying the accuracy of 
the solution. This dissipation operator is constructed in the following manner: 


Dz(Q) d i-^J t k 


(4.13) 
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r . 2 PiJ,k+Pi-lJ,k\ 

,,J ’ IP*+1 j,k + 2 J>» j,fc + Pi-1 j,fcl 

Typical values for the second and fourth difference damping constants are: 


(4.17) 


2 _ 1 4 1 

4 " 64 


(4.18) 


The complete dissipation operator ^ is constructed as the sum of the dissipation 
operators in each of the respective coordinate directions as: 


D ij,k - + ( D ')ij,k + ( D e)i,j,k 


(4.19) 


4.5 Implicit Residual Smoothing 

An effective technique for extending the stability limit of many explicit schemes is 
implicit residual smoothing. Residual smoothing attempts to accelerate the propaga- 
tion of changes in the dependent variables by filtering the residuals of the calculation 
(which may also be interpreted as the local time derivative of the computational so- 
lution) at each time step. By enhancing the transfer of information between grid 
points, calculation time steps much larger than the stability-limited values may be 
utilized. Residual smoothing was originally introduced by Lerat (see e.g. Hollanders, 
et al. [19]) for use with the Lax-Wendroff scheme and later applied to Runge-Kutta 
schemes by Jameson and Baker [20] as a technique to acclerate convergence for steady- 
state calculations. Jorgenson and Chima [21] applied this technique for unsteady flows 
and obtained a speedup of about 5 for their calculations without any change in the 
unsteady solution. 
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Since the time rate of change of the dependent variables dQ/dt is in essence 
controlled by a residual operator R(Q) = L(Q ) — D(Q), it would follow that any 
measure which accelerates the propagation of changes in the residual throughout 
the domain would ultimately enhance convergence. The implicit residual smoothing 
operator used in this study can be written as: 

(1 - e z 6 Z z){ 1 - «r^rr)(l - jjt = R iJ,k (4.20) 

where the differencing operator 6 is expressed as: 

8zzR iJ,k = R i+l,j,k ~ ^ R i,j,k + R i—lJ,k (4-21) 

A value of e Z z = €rr = = 2 is typically used. 

The reduction is applied along each coordinate direction separately as: 


R i,j,k = t 1 “ £zSzz ) lR iJ,k 

( 4 . 22 ) 

R i*j,k = C 1 “ €z 6 rr )~ lR iJ,k 

( 4 . 23 ) 

R i*j*k = ( 1_ ez 8 ee)~ lR **j,k 

( 4 . 24 ) 

n t>*** 

R i,j,k ~ R i,j,k 

( 4 . 25 ) 


where each of the first three steps above requires the inversion of a scalar tridiagonal 
matrix. The residual smoothing operator is applied at the first and third stage during 
the four-stage Runge-Kutta algorithm. The time-marching scheme then becomes: 

<?i = Q n - °i A{Q n ) 


Qi = Q"~ 
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Qz = Q n - <*3#(Q2) 

Q 4 = Q n - a A R(Q z ) 

Q n+1 = Q a (4.26) 

The implicit residual smoothing operator applied in this context can increase the 
effective time step by a factor of 3 or more, and can result in a 50% reduction in 
overall CPU time. 


4.6 Boundary Conditions 


Inflow and exit boundary conditions are applied numerically using characteristic 
theory. A one-dimensional isentropic system of equations is utilized to derive the 
following characteristic equations at an axial inflow/outflow boundary: 

dc~ , , . 

__ (t , z _ o )_ = 0 (4.27) 

dC+ , ,3C+ „ 

-^ + (»*+<>)— = 0 (4.28) 

where: 

2 a | 2a 

C — v z C + =v z + -1— (4.29) 

7-1 7-1 v ' 

In order to efficiently process boundary information in the numerical solution, phan- 
tom cells are located just outside the computational domain to permit the unmodified 
application of the interior point scheme at near boundary cells. Boundary condition 
information is effectively introduced into the solution by properly controlling the 
dependent variables in the phantom cells while permitting the application of the 
standard interior point scheme at near boundary cells. 
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For subsonic normal inflow, the upstream running invariant C~ is extrapolated 
to the inlet, and along with the equation of state, specified total pressure, total 
temperature, and flow angle, the flow variables at the boundary may be determined. 
A second boundary condition was available as an option which specifies the mass flow 
through the inlet, although this option is not recommended. At the exit, a static 
pressure is specified at the hub for internal flows, and at the outer boundary for 
external flows. The remaining pressures along the outflow boundary are calculated 
by integrating the radial momentum equation: 


dp PVff 
dr r 


( 4 . 30 ) 


In this case, the downstream running invariant C+ is used to update the phantom 
cells at the exit boundary. Far-field boundaries also use this characteristic technique 
based on whether the local flow normal to the boundary passes into or out of the 
domain. The solid surfaces (hub, cowl, airfoils) must satisfy flow tangency: 


V -n = 0 


( 4 . 31 ) 


In this case, we specify no convective flux through the boundary (an impermeable 
surface), and hence, only pressure is needed at the phantom cell. The pressure may 
be extrapolated, or updated using a variant of the normal momentum equation. In 
this case, extrapolation was found to be the most effective technique based on rapid 
convergence and adequate results. 
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4.7 Multiple-Block Coupling 


For the multiple-block C-grid scheme, the Euler solution is performed on a single 
grid block at a time. Special boundary conditions along block boundaries are therefore 
required to provide some transport of information between blocks. This transport is 
provided through a simple procedure which relies on the fact that the grid block 
boundaries have coincident grid points. Since the standard interior point scheme 
is performed at the first grid cell off the block boundary, all that is required is to 
determine the fluxes along the block boundaries themselves. Phantom points are 
provided along each side of a block boundary to accommodate the flux calculation at 
the block boundary cell face. The flow variables in the phantom cells are provided 
by interrogating the corresponding cell in the adjacent block when a particular value 
is required. After each stage of the Runge-Kutta integration for each block, the 
block boundary phantom cells are updated by utilizing the new dependent variable 
values from the adjacent blocks. This direct specification technique was compared 
to a characteristic-based method similar to the inlet and exit boundary condition 
routines, and was found to provide enhanced convergence behavior. 

Artificial damping is applied at the block boundaries by neglecting the fourth 
order derivative term. Implicit residual smoothing is applied at the block boundary 
by imposing a zero residual gradient (i.e. ( dR/dz ) = 0.0) condition at the boundary. 
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4.8 Solution Procedure 

The numerical solution is performed in an identical manner for both the H-grid 
and multiple-block C-grid Euler analyses. Assuming that the numerical grid and flow 
parameters are known, the time-marching procedure may begin from some set of ini- 
tial data. This initial data is initially specified as a uniform flow, or may be implied 
from a previous solution. The time-marching procedure is applied iteratively to up- 
date the flow variables as the solution proceeds. The solution is deemed converged 

O 

when the average residual R has been reduced by a factor of 10 

In the case of a counter-rotation calculation, the overall process of obtaining a 
converged solution is applied iteratively for each blade row as the source terms for the 
adjacent blade rows are updated. Thus, the time-marching loop is contained in an 
overall cycle which updates these source terms. The counter-rotation solutions were 
deemed converged when the differences in residuals between cycles were less than 

io -2 . 
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5. RESULTS 


Several numerical results from the 3D Euler analyses described in Chapter 3 are 
presented in the sections which follow. Initial results are given to verify the accuracy 
of the 3D Euler formulation for steady flow about unducted propfans. Calculations are 
presented for both the 8-bladed SR7 propfan geometry, and a 2-bladed SR7 propfan 
for which experimental airfoil surface static pressure distribution data were available. 
Several calculations for ducted propfans are presented. The numerics of the ducted 
propfan Euler analyses are first verified through comparisons with a two-dimensional 
numerical prediction of a circular arc cowl about a zero thickness nonrotating propeller 
aligned with the incoming flow. A second demonstration of the analyses is then 
performed for a fictitious ducted SR7 propfan. Predicted results are compared with 
numerical data obtained from a frequency domain panel method analysis for ducted 
propellers. A third set of calculations is given for a low speed ducted propeller. The 
predicted cowl surface pressure coefficient distribution is compared with experimental 
data for this geometry. Finally, the analyses are extensively compared with high 
speed experimental data for a 1.15 pressure ratio ducted fan configuration. Each case 
is discussed separately in the sections which follow. 
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5.1 SR7 8-Bladed Propfan 

In order to verify the numerics of the 3D Euler solver, several initial calculations 
were performed for steady flows about unducted propfans. A number of experimental 
studies have been performed for the unducted SR7 propfan geometry ranging from 
scaled wind tunnel tests to in-flight measurements. The SR7 propfan utilizes 8 blades 
with 41 degrees of sweep at the tip. A special hub contour is employed to eliminate 
flow choking at the hub. A description of the SR7 design parameters is given in 
Fig. 5.1. A particularly useful set of experimentally determined performance data 
for the SR7 was given by Stefko et al. [22]. In this study, a 2/9 scale model of the 
SR7 propfan was tested in a wind tunnel over a wide range of blade setting angles 
and operating conditions. A 70x35x15 (axial x radial x circumferential) numerical 
H-grid (as shown in Fig. 5.2) was generated to simulate this flow. Calculations were 
performed over a wide range of advance ratios and blade setting angles for flight Mach 
numbers of 0.7 and 0.8. All of the calculations for the unducted propfans described 
below were performed using the 3D H-grid Euler analysis operating in the unducted 
mode (no cowl boundary conditions are applied). 

A comparison of the predicted and experimental power coefficient map for the 
SR7 propfan at a flight Mach number of 0.8 is given in Fig. 5.3. During the investiga- 
tion of the SR7 geometry, it was immediately observed that the predicted results were 
extremely sensitive to the blade setting angle and required an accurate assessment 
of the deflected blade shape. Since no reasonable means was available to adjust the 
blade shape for each run, the blade angle was adjusted in a manner which permitted 
the predicted power coefficient to match the experimental power coefficient at the 
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design operating conditions. All other runs were made relative to this blade setting 
angle. The predicted power coefficient map illustrated in Fig. 5.3 is in fairly good 
agreement with the experimental curves for each blade setting angle. The comparison 
appears to worsen for larger blade setting angles as the power coefficient was consis- 
tently underpredicted by about 10%. This may be due to the inviscid flow assumption 
in the present analysis. 

A plot of the predicted propfan surface static/total pressure ratio contours for a 
Mach number of 0.8, 3/4 blade setting angle of 60.2 degrees, and an advance ratio of 
3.06 is given in Fig. 5.4. This contour plot clearly illustrates the spinner and blade 
root stagnation regions, and a shock which extends along the suction surface of the 
propfan blade. 

A similar series of calculations was performed at a flight Mach number of 0.7. 
In this case, an examination of the blade spanwise loading profiles was performed. A 
comparison of the predicted spanwise blade elemented power coefficient distributions 
with experimental data is given for several loading factors for M = 0.7 for the 8-bladed 
SR7 propfan in Fig. 5.5. The predicted power loading factors were somewhat different 
from the experimentedly measured values; however, it is obvious that the experimental 
data brackets the predicted spanwise distributions for each power loading value. 

5.2 SR7 2-Bladed Propfan Modane Tests 

A second, more extensive comparison of predicted results for steady unducted 
propfan flows was performed based on the test of the SR7 airfoil at the Modane wind 
tunnel test facility reported by Bushnell [23]. In this case, the propfan driver did 



Characteristic 


SR7 


Number of blades 

Tip sweep angle/ degrees 

Tip Speed (m/s) 

Power loading (kW/m**2) 

Activiy factor 

Integrated design lift coefficient 
Airfoils 

Ratio of nacelle maximum diameter to propeller diameter 
Cruise design Mach number 
Cruise design advance ratio 
Cruise design power coefficient 


8 

41 

243.8 

256.85 

227 

0.202 

NACA 16 and 65/CA 
0.35 
0.8 
3.06 
1.45 


Figure 5.1: SR7 propfan design characteristics 
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Figure 5.2: Axisymmetric plane projection of 70x35x15 SR7 grid 
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SR7 Power Coefficient Performance (M=0.8) 
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Figure 5.3: Comparison of predicted, experimental, and design power coefficient dis- 

tributions for 8-bladed SR7 propfan (M=0.8) 
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□ □ Experiment: Power Loading (Cp/J**3) 74.1% 

a a Experiment: Power Loading (Cp/J**3) 95.7% 

v v Experiment: Power Loading (Cp/J**3) 124.0% 

g o Experiment: Power Loading (Cp/J**3) 147.0% 


3D Euler: Power Loading (Cp/J**3) 66.0% 
3D Euler: Power Loading (Cp/J**3) 109.0% 



Figure 5.5: Comparison of predicted and experimental blade spanwise elemental 

power coefficient loading distributions for the 8-bladed SR7 propfan 
(M=0.7) 
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not have enough power to drive the 8-bladed configuration, so a 2-bladed version was 
tested instead. The airfoil surfaces were instrumented to permit the measurement of 
steady airfoil surface static pressures. The geometry and 120x45x25 H-grid used for 
these calculations are given in Fig. 5.6. This grid system utilized 49 points axially and 
21 points radially on the airfoil surface, which should be considered near minimum 
to accurately predict airfoil surface pressure distributions. 

Predictions were compared for a 3/4 radius blade setting angle of 30.58 degrees. 
The freestream Mach number was 0.2, with an advance ratio of 0.881. Predicted airfoil 
surface static pressure coefficient distributions are compared with experimental data 
at spanwise locations of 9.9%, 55.6%, and 94.4% (r/72) in Figs. 5. 7-5.9, respectively. 
The agreement between experiment and prediction is generally very good, with some 
small disagreement near the leading edge of each cross section. This discrepancy may 
be due to several factors. As in the 8-bladed SR7 calculations, the exact deflected 
airfoil shape and blade setting angle are unknown (experimental uncertainty is on the 
order of 1 degree), which introduces a significant source of error in the calculation. 
In addition, this particular flow condition gives rise to a leading edge vortex which 
bends across the suction surface of the blade. This feature is evident in both the 
experimental data and the calculation. 

In order to qualitatively compare the predicted airfoil surface pressure distribu- 
tions, a comparison of the predicted and experimental airfoil surface static pressure 
contours is given in Figs. 5.10-5.11 for the pressure (face) and suction (camber) sur- 
faces, respectively. The pressure surface shows little variation; however, the suction 
surface displays a noticable region of high pressure (presumably a result of the leading 
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edge vortex) which bends across the airfoil surface near the leading edge from roughly 
50%-100% span. This feature is also evident in the experimental static pressure con- 
tours, and illustrates the remarkable similarity between prediction and experiment. 


5.3 Symmetric Circular Arc Cowl Test Case 

In order to verify the accuracy of the Euler formulation for a ducted propfan 
geometry, a numerical test case was initiated based on a symmetric double circular 
arc airfoil cowl about a zero thickness, nonrotating propeller with a large hub diame- 
ter. A single passage of this geometry is given in Fig. 5.12. The propeller blades are 
aligned with the incoming flow and are therefore “invisible” to the inviscid numerical 
solution. These aspects, combined with the large hub diameter and large number of 
propeller airfoils (320) essentially reduce the aerodynamics to two dimensions. There- 
fore, the 3D Euler results were compared with an extensively verified 2D Euler code 
based on the hopscotch [24] time-marching technique for a similar two-dimensional 
geometry. A comparison with numerical, rather than experimental, data should be 
more enlightening in this instance since discrepancies between experimental data and 
calculations may not necessarily be due to boundary condition errors (i.e., viscous 
effects, grid induced errors, unsteadiness, etc.). 

The 3D H-grid calculations were performed on a 91x41x5 grid, while the 2D 
calculations were performed on a 91x21 grid. A comparison of the 2D grid and 
the 3D grid axisymmetric projection is given in Fig. 5.13. Since the 2D problem is 
perfectly symmetric, only half of the flowfield was computed. Each grid utilized 31 
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Figure 5.6: 2-bladed SR7 propfan geometry and 120x45x25 grid for Modane test 
comparison 
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2-Bladed SR7 Modane Test Comparison (M=0.2, 9.9% span) 



Axial Chord (x/Cx) 


Figure 5.7: Comparison of predicted and experimental airfoil surface static pressure 

coefficient distributions for 2-bladed SR7 propfan Modane test (9.9% 
span) 
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Figure 5.11: Comparison of predicted and experimental airfoil suction surface 

static/total pressure ratio contours for 2-bladed SR7 propfan 
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points to define the cowl upper &nd lower surfaces. 

Calculations were compared for inflow Mach numbers of 0.5 and 0.7. A com- 
parison of the predicted cowl surface grid line Mach number distributions is given in 
Fig. 5.14. The 3D Euler results for both the inner and outer cowl surfaces display 
excellent agreement with the 2D Euler code results. A similar comparison is given in 
Fig. 5.15 for the 0.7 Mach number case. For this case, a strong symmetric shock wave 
pattern forms on the cowl surface, affording an examination of the shock- capturing 
characteristics of the cowl boundary algorithm. In this case, small discrepancies exist 
between the 2D solution, and the upper and lower surfaces of the 3D solution, al- 
though the agreement is still generally very good. The rapid decrease in Mach number 
across the shock wave is evident in each calculation. The small discrepancies between 
the solutions, and between the upper and lower surfaces of the 3D Euler solution, are 
easily explained by the differences in the mesh spacing near the cowl. 

A similar comparison was performed for the multiple-block C-grid Euler solver. 
In this case, the five-block grid arrangement shown in Fig. 5.16 was generated for the 
numerical predictions. Again, the 91x21 2D grid is also displayed for comparison. 
The total number of grid points is roughly similar to the H-grid calculations with 
individual block sizes of (65x13x5), (16x13x5), (81x7x5), (16x13x5), and (65x14x5). 
A comparison of the predicted cowl surface grid line Mach number distribution horn 
the 3D multiple-block C-grid Euler analysis and the 2D solution for an inflow Mach 
number of 0.5 is given in Fig. 5.17. Again, good agreement is observed between 
the 3D and 2D results. A second comparison of results is given for an inflow Mach 
number of 0.7 in Fig. 5.18. Again, the predicted Mach number distributions are in 
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good agreement, and the effect of the shock which results from the higher inflow 
Mach number is clearly displayed. Static pressure contour plots for the axisymmetric 
plane projection of the 3D H-grid and multiple-block C-grid Euler results for the 0.7 
inflow Mach number case are given in Fig. 5.19. The individual block boundaries are 
outlined in this case to demonstrate the smooth transition of contour lines across the 
block interfaces. The H-grid and multiple-block C-grid results are very similar in spite 
of the differences in the solution technique. The asymmetry of the contour lines is a 
result of the differences in boundary conditions at the hub and outer boundaries. This 
figure serves to validate the accuracy of the block coupling procedure and indicates 
that no nonphysical behavior is introduced by this procedure. 

5.4 Ducted SR7 Demonstration 

An interesting comparison of results was based on the calculation of the flowfield 
about a fictitious ducted 8-bladed SR7 propfan. In spite of the fact that there are 
no experimental data for this configuration, these calculations were useful for pre- 
liminary debugging of the analysis for ducted configurations, and also provided an 
opportunity to compare the predicted Euler results with an existing frequency domain 
panel method also developed for ducted propfan flows. 

The geometry was based on the original 8-bladed SR7 geometry given in Section 
4.1 with a thin profile duct about the propeller. A relatively large tip gap (2% of the 
duct diameter) was imposed to permit an examination of the tip flow. 

The geometry and grid surfaces are shown in Fig. 5.20. The axisymmetric pro- 
jection of the 69x45x15 H-grid generated for this case is illustrated in Fig. 5.21. The 
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Figure 5.12: Symmetric double circular arc cowl ducted propeller test case geometry 
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2D 



3D 



Figure 5.13: Comparison of 91x21 2D grid and axisymmetric plane projection of 

91x42x5 3D H-grid for the symmetric double circular arc cowl test case 





Circular Arc Cowl Test Case (M=0.7) 



Figure 5.15: Comparison of 2D and 3D H-grid predicted cowl surface grid line Mach 

number distributions (M=0.7) 
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Figure 5.16: Comparison of 91x21 2D grid and axisymmetric plane projection of 
multiple-block C-grid for the symmetric double circular arc cowl test 


case 
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Circular Arc Cowl Test Case (M=0.5) 



Figure 5.17: 


Comparison of 2D and 3D multiple-block C-grid predicted cowl surface 
grid line Mach number distributions (M=0.5) 




Figure 5.19: Comparison of the predicted static pressure contours for the axisym- 

metric plane projection of the 3D H-grid (top) and multiple-block C- 
grid (bottom) Euler analyses of the symmetric double circular arc cowl 
test case (M=0.7) 
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case tested was for a true 3/4 blade setting angle of 60.2 degrees, flight Mach num- 
ber of 0.7 , and an advance ratio of 3.06. Predicted results from the 3D Euler code 
were compared with predicted results from a frequency domain panel code developed 
for ducted propellers by Williams et al. [25]. A comparison of the predicted airfoil 
spanwise elemental power coefficient loading distributions is given in Fig. 5.22 for 
this case. A similar comparison is given for the predicted airfoil spanwise elemental 
thrust coefficient loading distributions in Fig. 5.23. The predicted power and thrust 
coefficient loading distributions for the unducted case are also given in Figs. 5.22-5.23 
for comparison. Unfortunately, the panel code does not accurately model the hub 
and cowl contour; however, the loading distributions are largely similar, although the 
Euler code predicts a somewhat larger aerodynamic loading near the blade tip region. 
The predicted propfan surface static pressure contours for the 3D Euler results are 
pictured in Fig. 5.24. The effect of adding the duct to the normally unducted blade 
is immediately evident in the predicted results. The ducted blade becomes highly 
loaded at the tip, and generates a rather strong shock which extends to the inner 
duct surface. The spanwise loading indicates that the addition of the duct results in 
a 30% increase in thrust and a 25% increase in power. 

A second 3D solution was performed for this geometry using the multiple-block 
C-grid Euler analysis. In this case, the block grid pictured in Fig. 5.25 was utilized. 
The individual grid blocks contained (69x21x15), (11x9x15), (103x5x15), (12x9x15), 
and (69x18x15) points, respectively. Again, the predicted results were compared to 
the frequency-domain panel method of Williams et al. [25]. A comparison of the 
predicted airfoil spanwise elemental power coefficient loading distributions is given in 
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Fig. 5.26 for this case. A similar comparison is given for the predicted airfoil spanwise 
elemental thrust coefficient loading distributions in Fig. 5.27. The results exhibited 
by the multiple-block C-grid Euler analysis are largely similar to the H-grid Euler 
analysis results. 


5.5 Low Speed Ducted Propeller 

The third ducted configuration tested was based on the low speed ducted pro- 
peller originally tested by Kruger [26]. Kruger’s geometry was based on an axisym- 
metric elliptical centerbody with an 8-bladed NACA 23009-23012 (hub to tip) section 
propeller. Data were taken for the static pressure distributions on the centerbody and 
cowl surfaces. 

A 160x41x15 grid was generated for this geometry, as shown in Fig. 5.28. The 
axisymmetric projection of this grid is given in Fig. 5.29. The case presented corre- 
sponds to Kruger’s cowl #1, propeller #1, and nacelle #1. The calculation was based 
on an advance ratio of 0.95 and a flight Mach number of 0.4275. This experimental 
arrangement was particularly convenient for numerical simulation since the duct was 
supported separately from the centerbody, and thus mounting struts which might 
cause aerodynamic interference were not required. 

A comparison of the H-grid predicted and experimental cowl surface static pres- 
sure coefficient distributions is given in Fig. 5.30. The somewhat jagged appearance 
of the computational results is due to a lack of smoothness in the coordinates used to 
define Kruger’s cowl geometry (the cowl coordinates were obtained by manually digi- 
tizing a drawing of the cowl cross section). In spite of the crudeness of the geometric 
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Figure 5.20: Ducted 8-bladed SR7 geometry and 69x45x15 H-grid 
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Figure 5.21: Axisymmetric plane projection of 69x45x15 H-grid for ducted SR7 

propfan 
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Ducted SR7 Spanwise Load Analysis M=0.7 J=3.06 



Figure 5.22: Comparison of H-grid predicted elemental blade power coefficient distri- 

bution with frequency domain panel method results for ducted 8-bladed 
SR7 
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Ducted SR7 Spanwise Load Analysis M=0.7 J=3.06 



Figure 5.23: 


Comparison of H-grid predicted blade elemental thrust coefficient dis- 
tributions with frequency domain panel method results for ducted 8- 


bladed SR7 



Static Pressure Ratio 
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Figure 5.24: Predicted H-grid propfan surface static/total pressure ratio contours 

for ducted SR7 
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Figure 5.25: Axisymmetric plane projection of multiple-block C-grid for ducted SR7 

propfan 
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Ducted SR7 Spanwise Load Analysis M=0.7 J=3.06 



Figure 5.26: Comparison of multiple-block C-grid predicted elemental blade power 
coefficient distribution with frequency domain panel method results for 
ducted 8-bladed SR7 
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Ducted SR7 Spanwise Load Analysis M=0.7 J=3.06 
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Comparison of multiple-block C-grid predicted blade elemental thrust 
coefficient distributions with frequency domain panel method results 
for ducted 8-bladed SR7 
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Figure 5.28: Geometry and 127x48x15 grid for Kruger low speed ducted propeller 
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Figure 5.29: Axisymmetric plane projection of 127x48x15 H-grid for Kruger low 

speed ducted propeller 
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model, good agreement was obtained with the experimental data. 

A second calculation was performed for this case using the 3D multiple-block C- 
gnd Euler analysis. The axisymmetric plane projection of the calculation grid is given 
in Fig. 5.31. The respective grid block sizes for this calculations were (150x21x15), 
(65x13x15), (147x7x15), (19x13x15), and (150x15x15). Again, a comparison of the 
multiple-block C-grid predicted and experimental cowl surface local static pressure 
coefficient distributions is given in Fig. 5.32. A similar level of agreement was obtained 
from the multiple-block C-grid calculation, in spite of the geometric uncertainties. 


5.6 NASA 1.15 Pressure Ratio Fan 


Final calculations were performed for a 1.15 pressure ratio fan stage originally 
tested at NASA for a wide range of flows [27]-[30]. A description of the geometry is 
given in Fig. 5.33. This geometry is representative of a 25:1 bypass ratio turbofan en- 
gine fan stage. Both single-rotation and counter-rotation calculations were performed 
for this geometry. For each calculation, it was not possible to exactly match the ex- 
perimentally measured mass flow through the fan due to the expense of iterating on 
the three-dimensional solution; therefore, predictions are based on the flight Mach 
number and estimated fan rotational speed alone. A three-dimensional view of the 
geometry and grid system overlay are given in Fig. 5.34. 

H-grid predictions were based on the 225x52x15 H-grid pictured in the axisym- 
metric plane projection given in Fig. 5.35. The blades were represented by 23 points 
in the axial direction and 21 points in the radial direction. Both blade rows are rep- 
resented in the axisymmetric projection of the grid. Single-rotation calculations were 
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Kruger Ducted Propeller Test (M=0.4275, J=0.95) 



Figure 5.30: Comparison of H-grid predicted and experimental cowl surface static 

pressure coefficient distributions for Kruger low "speed propeller 
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Figure 5.31: Axisymmetric plane projection of multiple-block C-grid for Kruger low 

speed ducted propeller 
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Kruger Ducted Propeller Test (M=0.4275, J=0.95) 



Figure 5.32: Comparison of multiple-block C-grid predicted and experimental cowl 

surface static pressure coefficient distributions for Kruger low speed 
propeller 
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Figure 5.33: NASA 1.15 pressure ratio fan stage geometry (dimensions in cm) 
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Figure 5.34: NASA 1.15 pressure ratio fan stage and 3D grid overlay 
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performed by neglecting the thickness of the deswirl vanes. 

For the initial calculations, only the single-rotation rotor was modelled. A com- 
parison of predicted and experimental cowl leading edge surface static pressure co- 
efficient distributions for a flight Mach number of 0.75 and an advance ratio ( J) of 
2.86 for the rotor alone is given in Fig. 5.36. The data are presented in an unwrapped 
form to expand the critical region of interest at the cowl leading edge. The agree- 
ment between prediction and experiment is outstanding. Each of the features of 
the experimental pressure distribution is well defined by the predicted results. Any 
discrepancies between experiment and prediction are thought to be due to small dif- 
ferences between the predicted and experimentally measured mass flow rate through 
the fan or grid inadequacies. A similar comparison of results is given in Fig. 5.37 for a 
flight Mach number of 0.85 and an advance ratio of 3.23. In this case, good agreement 
is observed over much of the cowl surface, except for a rather large discrepancy on 
the outer cowl surface downstream of the leading edge region. Experimental evidence 
indicates that a flow separation occurs for the outer cowl surface in this case, and the 
disagreement is a result of the inviscid flow assumption in the present analysis. This 
would suggest that a viscous flow analysis will ultimately be required for high speed 
ducted fan design and off-design analysis. 

In order to assess the impact of the deswirl vanes on the aerodynamic solution, the 
0.75 Mach number solution was repeated for a counter-rotation configuration using 
the average-passage approach. This demonstration was performed under a no-cost 
contract amendment to verify the average-passage formulation in the ducted propfan 
Euler solver. The solution procedure involves iteratively updating representative body 
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Figure 5.35: 


Axisymmetric plane projection of NASA 1.15 pressure ratio fan stage 
225x52x15 H-grid 
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1.15 Pressure Ratio Fan Cowl Pressure Coefficient (M=0.75) 



Area Ratio (A/A(MAX)-1 .0) 


Figure 5.36: Comparison of predicted and experimental cowl leading edge static 

pressure coefficient distributions for NASA 1.15 pressure ratio fan (H- 
grid, rotor only, M=0.75) 
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forces for the adjacent blade in each individual blade row solution. The introduction 
of the axisymmetric cowl has no bearing on this solution procedure. 

Again, the grid pictured in Fig. 5.35 was used although now the thickness of the 
deswirl vanes is included. A comparison of predicted and experimental data for the 
cowl leading edge surface static pressure coefficient distributions, nozzle boattail sur- 
face static pressure coefficient distributions, and nozzle inlet radial total pressure ratio 
distributions is given in Fig. 5.38. Again good agreement is obtained for each of these 
stations, demonstrating the capability of accurately simulating the aerodynamics of 
a complete ducted fan engine configuration. 

An illustration of the predicted fan surface static pressure ratio contours for this 
calculation is given in Fig. 5.39. The extreme flow gradients at the cowl leading edge 
and the aerodynamic interaction between the blade tips and cowl inner surface are 
well displayed in the predicted contours. 

The single-rotation fan geometry flowfield was also calculated using the multiple- 
block C-grid Euler solver. Predictions were based on the multiple-block grid pictured 
in Fig. 5.40. Grid block sizes were generated to be compatible with the total number 
of points used in the H-grid calculation, with block sizes of (172x25x15), (68x17x15), 
(193x9x15), (17x17x15), and (172x21x15), respectively. A comparison of the pre- 
dicted and experimental cowl leading edge surface static pressure ratios is given in 
Fig. 5.41 for a Mach number of 0.75 and an advance ratio of 2.86. Again, excellent 
agreement between experiment and calculation is observed. The embedded C-grid 
block permits a large number of grid points in the region of the cowl leading edge 
without incurring an excessive number of points in the remainder of the flowfield. The 
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Figure 5.38: Comparison of predicted and experimental pressure data for NASA 1.15 

pressure ratio fan stage (H-grid, full stage, M=0.75) 


Static Pressure Ratio 


ORIGINAL PAGE 
COLOR PHOTOGRAPH 


91 




92 


high shear associated with the H-grid in this region is also eliminated. The solution 
was found to have a smooth continuous transition across block boundaries, and a 
relatively rapid convergence (1500 iterations) was achieved. 

A second calculation was performed for the C-grid algorithm at a Mach number of 
0.85. A comparison of the predicted and experimented cowl leading edge surface static 
pressure coefficient distributions is given in Fig. 5.42. Again, excellent agreement 
is observed over most of the cowl surface, except near the downstream portion of 
the cowl outer surface. This is again thought to be due to an experimental flow 
separation, although the deviation is less pronounced in this case. The C-grid results 
do not differ significantly from the H-grid results for either of the cases examined, 
although a detailed grid resolution study has not been performed. In addition, the 
H-grid calculation was performed with a relatively large number of grid points about 
the cowl leading edge to enhance the accuracy in this critical region. As the number 
of grid points decreases, the C-grid analysis will likely be superior. 
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Figure 5.40: Multiple-block C-grid for NASA 1.15 pressure ratio fan stage (rotor 

only) 
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0. CONCLUSIONS 


A three-dimensional Euler analysis and grid generation scheme have been de- 
veloped for the numerical analysis of ducted propfan flowfields. Both H-grid and 
multiple-block C-grid solution structures have been examined. The analysis has been 
verified through comparisons with experimental results from both a low speed ducted 
propeller and a high speed 1.15 pressure ratio fan stage for both single-rotation and 
counter-rotation configurations. At least one case was observed which indicated that 
a viscous flow analysis might improve the accuracy of the predictions. A comparison 
of the predicted results for both the H-grid and multiple-block C-grid analysis did 
not display significant differences, although as the number of grid points decreases, 
the C-grid analysis will likely be better. The capability of accurately simulating the 
aerodynamics about a complete ducted fan propulsor system has been demonstrated. 

Clearly the largest obstacle in accurately predicting propfan flowfields of any 
kind is in the specification of the deflected airfoil shape. Calculations for the SR7 
clearly displayed how even small errors in the blade setting angle could drastically 
change the predicted results. A second obstacle encountered in this study was a lack 
of experimental data for ducted fans which prevented a lengthy verification study for 
a range of geometries and flow conditions. 
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APPENDIX A. GRID GENERATOR OPERATING INSTRUCTIONS 


A.l Introduction 

Grid generation for both the H-grid and C-grid Euler analyses is performed by a 
single program with identical input formats. The grid generation source program is 
written in FORTRAN 77, and has been used successfully on Cray UNICOS and IBM 
VM/CMS mainframe computer systems, as well as Silicon Graphics 4D workstations 
using a UNIX operating system. Array dimensions in the programs lire specified by 
a PARAMETER statement in the MAIN routine. The PARAMETER statement 
appears as: 

PARAMETER( IJX=100, KX=15, NX=2) 
and appears in every subroutine. Another important parameter, PRECIS, is based 
on the precision of the host computer and is determined internally. 

The individual variables in the PARAMETER statement are defined as: 

IJX Used to determine the maximum number of points in the axial direction 
KX The maximum number of points in the circumferential direction 
NX The maximum number of blade rows 

Array storage is determined through a variable IXT which is defined as IXT = 
(2 NX + 1)1 JX. The parameter IXT then determines the maximum number of 



104 


points in the axial direction. Several rules for storage must be observed for proper 
operation of the code. Storage is assigned under the assumption that IJX > KX. 
The parameter NX must be greater than or equal to the number of blade rows. For 
multiple-block C-grid generation, KX must be greater than or equal to JMAX + 1 
(JMAX is an input parameter described below). The first step in any grid generation 
sequence should be to estimate the size of the grid to be generated and to make sure 
that adequate array storage is provided for the anticipated grid. Unusual results can 
occur when an array size is exceeded, and many unexplained errors can be traced to 
inadequate array dimensions. 

The grid generation code can be used for either ducted or unducted, single- or 
counter-rotation propfans systems in the H-grid mode, and ducted single-rotation 
propellers in the C-grid mode. An option is provided to arbitrarily adjust the 3/4 
radius blade setting angle without modifying the blade geometry parameters. This 
assumes that the pitch change axis is perpendicular to the axis of rotation, and only a 
solid body rotation can occur (no deflections). This procedure is somewhat limited by 
the airfoil cross section coordinate distribution due to the requirement of maintaining 
the index of the true leading and trailing edge points. 

The program produces plotting output suitable for viewing using the PLOT3D 
graphics software package developed at the NASA Ames Research Center. This soft- 
ware permits workstation grid generation and processing before invoking the flow 
solver. PLOT3D output is written to FORTRAN unit 4 for H-grid generation, while 
the C-grid version writes individual block grid information to units 41-45 for the 
five-block architecture. Additional PLOT3D outputs can be generated for debugging 
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purposes by small changes in the input. Some PostScript output is also available on 
FORTRAN unit 15, although this is primarily for analyzing the quality of the cowl 
C-grid in the multiblock grid generation. 

The input/output files used by the grid generation codes are described in the 
first section below. This section is then followed by a more general section which 
describes the function of the individual subroutines in the computer program. 

Geometric data may be input in either dimensional, or nondimensional form. 
However, each component must be defined in a consistent (either dimensional or 
nondimensional) manner. The final grid coordinates will be nondimensionalized by 
the maximum propeller diameter. 

A. 2 Input/Output File Description 

An example input data file is shown in Fig. A.l. A brief description of the 
variables used in the data file is given below. 

VARIABLE DESCRIPTION 


TITLE An 80 character title for the mesh generation. 

NBLROW Number of blade rows: 

for C-grids, this must be 1, 

for H-grids, this must be less than or equal to PARAMETER NX. 
IGEOM Internal/external geometry parameter: 

if = 1, internal flow (e.g. compressor, turbine), 
if = 2, external flow (e.g. propfan, ducted propfan). 
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+TITLE 

GENERIC DUCTED PROPFAN 

+NBLROW — ++IGEOM ++ITHETA — ++IWRITE — ++NCOWL ++JMAX ++IORTH ++IDEBUG — 

1 2 2 2 1 5 1 

+NINPLS — ++NEXPLS — ++NBLPTZ — ++NBLPTR — ++NBLPTT — ++CITER ++ISOLVE — ++ 


1 

1 

15 

15 

7 

50 

0 


+ZINLET — ++ZEXIT ++CFACLE — ++CFACTE — ++DFACLE — ++DFACTE — ++- 

+ + - 


-0.50000 

1.09000 

0.02 

0.02 

1.0 

1.0 



+RATIN ++RATEX ++RATBB ++RATBLZ — ++RATBLR — ++RATBLT — ++RATTOB — ++- 


1.32000 

1.32000 

1.20000 

1.200000 

1.1000 

1.05000 

1.36000 



+NHUB ++NT IP ++ - 

+ + - 

++- 

++- 

+ + - 

+ + - 


53 

2 







+ZHUB,RHUB+- 

+ + - 





++- 

+ + - 


-80.00000 - 

■20.00000 

-0.45000 

-0.40000 

-0.37000 

-0.35000 

-0.33600 

-0.330C 

-0.32685 

-0.32336 

-0.31773 

-0.31097 

-0.30302 

-0.29383 

-0.28285 

-0.257$ 

-0.22232 

-0.18130 

-0.13606 

-0.09774 

-0.07757 

-0.07042 

-0.06853 

-0.061$ 

-0.05049 

-0.03837 

-0.02810 

-0.01930 

-0.01181 

-0.00522 

0.00080 

0.006$ 

0.01215 

0.01808 

0.02469 

0.03237 

0.04148 

0.05238 

0.06434 

0.073] 

0.07788 

0.07965 

0.09454 

0.11764 

0.15665 

0.21074 

0.27769 

0.354$ 

0.43958 

0.58799 

0.75179 

0.92019 

1.09000 




0.00200 

0.00200 

0.00200 

0.00200 

0.00260 

0.00400 

0.00600 

0.008C 

0.01000 

0.01300 

0.02177 

0.03076 

0.03877 

0.04635 

0.05436 

0.068] 

0.08364 

0.09408 

0.09971 

0.10164 

0.10204 

0.10239 

0.10251 

0.103C 

0.10421 

0.10611 

0.10819 

0.11031 

0.11257 

0.11463 

0.11653 

0.118$ 

0.12065 

0.12307 

0.12572 

0.12874 

0.13223 

0.13630 

0.14061 

0.143( 

0.14530 

0.14590 

0.15081 

0.15819 

0.16767 

0.17458 

0.17754 

0.174$ 

0.16806 

0.15661 

0.15465 

0.15465 

0.15465 




+ZTIP , RTIP+- 


++- 

+ +- 

+ +- 

+ +- 

++- 



- 10.0000 10.000 
1.7000 1.700 

******************************************************************************’ 
+NSHR ++NSHRIN — ++NSHROU — ++NPBCAB — ++ ++ ++ ++ 


75 

7 

7 

3 





+ ZSHR,RSHR+ 

+ + 

+ + 

+ + 

++ 

+ + 

+ + 


0.250000 

0.249464 

0.247860 

0.245194 

0.241480 

0.236734 

0.230970 

0.2242: 

0.216504 

0.207864 

0.198340 

0.187960 

0.176774 

-0.164834 

-0.176774 

-0.18791 

-0.198340 

-0.207864 

-0.216504 

-0.224220 

-0.230970 

-0.236734 

-0.241480 

-0.2451! 

-0.247860 

-0.249464 

-0.250000 

-0.249464 

-0.247860 

-0.245194 

-0.241480 

-0.2367$ 

-0.230970 

-0.224220 

-0.216504 

-0.207864 

-0.198340 

-0.187960 

-0.176774 

-0.1648; 

-0.152190 

-0.138894 

-0.125000 

-0.110570 

-0.095670 

-0.080360 

-0.064706 

-0.0487' 

-0.032630 

-0.016350 

0.000000 

0.016350 

0.032630 

0.048774 

0.064706 

0.08031 

0.095670 

0.110570 

0.125000 

0.138894 

0.152190 

0.164834 

0.176774 

0 .1879* 

0.198340 

0.207864 

0.216504 

0.224220 

0.230970 

0.236734 

0.241480 

0.2451! 

0.247860 

0.249464 

0.250000 






0.410000 

0.410010 

0.410040 

0.410067 

0.410083 

0.410075 

0.410047 

0.41001 

0.409958 

0.409910 

0.409857 

0.409797 

0.409725 

0.404410 

0.404362 

0.40431 

0.404480 

0.404638 

0.404887 

0.405247 

0.405617 

0.406077 

0.406807 

0.4074- 

0.408125 

0.408895 

0.410000 

0.410750 

0.411665 

0.413025 

0.414677 

0.4161! 

0.417597 

0.419175 

0.420660 

0.422092 

0.423482 

0.424795 

0.426020 

0.4271! 


Figure A.l: Sample input data file for grid generation 
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0.428175 

0.430775 

0.423212 

0.414738 

0.410162 


0.429077 

0.430250 

0.422045 

0.413900 

0.410040 


0.429845 

0.429552 

0.420887 

0.413117 

0.410000 


0.430467 

0.428700 

0.419752 

0.412395 


0.430922 

0.427725 

0.418655 

0.411742 


0.431187 

0.426660 

0.417597 

0.411178 


0.431250 

0.425538 

0.416587 

0.410720 


0.431113 

0.424382 

0.415635 

0.410382 


********************************************* # ^^ # ^^^^ ## ^ # ^^^ 
+NBLD ++NBLCRC — ++NPPRC ++IPCH ++ZPCA ++THPCA ++BETA ++ + 

® 8 21 1 0.0 0.0 602 
+ ZBLA, RBLA, THBLA, TTBLA ++ ++ ++ ++ ---++ + 


- 0.07042 

- 0.00878 

0.05004 

- 0.07840 

- 0.01513 

0.04530 

- 0.08927 

- 0.02326 

0.04005 

- 0.10222 

- 0.03235 

0.03526 

- 0.11092 

- 0.03762 

0.03396 

- 0.10567 

- 0.03281 

0.03883 

- 0.08982 

- 0.02145 

0.04622 

- 0.06358 

- 0.00411 

0.05509 

0.10239 

0.11351 

0.13544 

0.11840 

0.12754 

0.14877 

0.14224 

0.14897 

0.16910 

0.17874 

0.18345 

0.20169 

0.22445 

0.22805 

0.24409 

0.27612 

0.27945 


- 0.06259 
- 0.00130 
0.05742 
- 0.07040 
- 0.00745 
0.05285 
- 0.08096 
- 0.01519 
0.04792 
- 0.09346 
- 0.02372 
0.04363 
- 0.10175 
- 0.02853 
0.04283 
- 0.09655 
- 0.02376 
0.04775 
- 0.08126 
- 0.01295 
0.05468 
- 0.05613 
0.00330 
0.06249 
0.10294 
0.11588 
0 . 13813 
0.11876 
0.12983 
0.15145 
0.14233 
0.15095 
0.17176 
0.17901 
0.18499 
0.20427 
0.22475 
0.22934 
0.24644 
0.27636 
0.28059 


- 0.05477 

0.00619 

0.06481 

- 0.06241 

0.00021 

0.06041 

- 0.07266 

- 0.00719 

0.05581 

- 0.08470 

- 0.01517 

0.05203 

- 0.09257 

- 0.01948 

0.05172 

- 0.08743 

- 0.01475 

0.05668 

- 0.07270 

- 0.00447 

0.06314 

- 0.04869 

0.01070 

0.06990 

0.10369 

0.11825 

0.14078 

0.11930 

0.13216 

0.15409 

0.14263 

0.15319 

0.17438 

0.17911 

0.18689 

0.20679 

0.22492 

0.23084 

0.24874 

0.27645 

0.28199 


- 0.04699 

0.01346 

0.07222 

- 0.05443 

0.00786 

0.06800 

- 0.06436 

0.00081 

0.06372 

- 0.07594 

- 0.00667 

0.06043 

- 0.08339 

- 0.01050 

0.06062 

- 0.07831 

- 0.00577 

0.06562 

- 0.06414 

0.00401 

0.07160 

- 0.04125 

0.01810 

0.07731 

0.10469 

0.12120 

0.14337 

0.12005 

0.13456 

0.15666 

0.14312 

0.15546 

0.17693 

0.17927 

0.18902 

0.20926 

0.22504 

0.23273 

0.25098 

0.27653 

0.28357 


- 0.03924 

0.02074 

0.07965 

- 0.04648 

0.01530 

0.07560 

- 0.05608 

0.00875 

0.07164 

- 0.06718 

0.00183 

0.06886 

- 0.07421 

- 0.00152 

0.06953 

- 0.06919 

0.00321 

0.07457 

- 0.05558 

0.01246 

0.08007 

- 0.03381 

0.02549 

0.08472 

0.10595 

0.12414 

0.14590 

0.12104 

0.13753 

0.15919 

0.14379 

0.15789 

0.17943 

0.17967 

0.19118 

0.21167 

0.22525 

0.23467 

0.25317 

0.27679 

0.28518 


- 0.03153 

0.02804 

- 0.03857 

0.02277 

- 0.04782 

0.01654 

- 0.05844 

0.01025 

- 0.06504 

0.00741 

- 0.06008 

0.01214 

- 0.04704 

0.02089 

- 0.02638 

0.03289 

0.10746 

0.12704 

0.12229 

0.14041 

0.14471 

0.16079 

0.18025 

0.19357 

0.22562 

0.23678 

0.27718 

0.28704 


- 0.02386 

0.03535 

- 0.03070 

0.03026 

- 0.03959 

0.02436 

- 0.04971 

0.01855 

- 0.05588 

0.01626 

- 0.05097 

0.02102 

- 0.03850 

0.02933 

- 0.01895 

0.04029 

0.10915 

0.12989 

0.12378 

0.14325 

0.14589 

0.16361 

0.18105 

0.19638 

0.22620 

0.23921 

0.27773 

0.28914 


- 0.01628 

0.04269 

- 0.02287 

0.03777 

- 0.03140 

0.03219 

- 0.04101 

0.02690 

- 0.04674 

0.02510 

- 0.04189 

0.02992 

- 0.02997 

0.03777 

- 0.01153 

0 . 04769 . 

0.11119 

0.13269 

0.12550 

0.14603 

0.14732 

0.16639 

0.18213 

0.19906 

0.22701 

0.24169 

0.27849 

0.29112 


Figure A.l: Sample input data file for grid generation (continued) 
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0.29306 

0.29497 

0.29683 

0.29863 

0.30040 




0.32968 

0.32972 

0.32988 

0.33011 

0.33045 

0.33092 

0.33153 

0 . 3322 < 

0.33312 

0.33418 

0.33533 

0.33651 

0.33789 

0.33939 

0.34084 

0.34221 

0.34365 

0.34500 

0.34632 

0.34761 

0.34886 




0.38143 

0.38163 

0.38192 

0.38227 

0.38269 

0.38317 

0.38374 

0 . 3844 : 

0.38513 

0.38585 

0.38668 

0.38761 

0.38851 

0.38940 

0.39028 

0 . 3911 : 

0.39196 

0.39278 

0.39357 

0.39436 

0.39512 




- 0.22623 

- 0.11448 

- 0.06580 

- 0.00593 

0.03783 

0.07747 

0.11208 

0.13841 

0.15770 

0.17116 

0.17992 

0.18148 

0.18524 

0.18762 

0.18808 

0 . 1867 : 

0.18416 

0.18144 

0.17914 

0.17788 

0.17801 




- 0.22180 

- 0.14123 

- 0.10573 

- 0.06270 

- 0.02902 

0.00160 

0.03009 

0 . 0550 : 

0.07632 

0.09372 

0.10771 

0.11843 

0.12448 

0.12892 

0.13273 

0 . 1371 ! 

0.14188 

0.14653 

0.15096 

0.15519 

0.15910 




- 0.21743 

- 0.16069 

- 0.13115 

- 0.09805 

- 0.07011 

- 0.04451 

- 0.02039 

0 . 0015 ! 

0.02138 

0.03888 

0.05401 

0.06720 

0.07870 

0.08821 

0.09698 

0.10501 

0.11278 

0.12067 

0.12867 

0.13656 

0.14430 




- 0.22009 

- 0.17785 

- 0.15193 

- 0.12440 

- 0.09990 

- 0.07691 

- 0.05491 

- 0 . 0343 ( 

- 0.01496 

0.00299 

0.01943 

0.03450 

0.04855 

0.06161 

0.07367 

0 . 0853 : 

0.09651 

0.10748 

0.11820 

0.12866 

0.13889 




- 0.22349 

- 0.18910 

- 0.16466 

- 0.13966 

- 0.11662 

- 0.09455 

- 0.07318 

- 0 . 0527 : 

- 0.03311 

- 0.01442 

0.00338 

0.02025 

0.03646 

0.05204 

0.06701 

0 . 0815 ! 

0.09562 

0.10941 

0.12280 

0.13576 

0.14805 




- 0.19509 

- 0.16591 

- 0.14364 

- 0.12121 

- 0.09998 

- 0.07939 

- 0.05926 

- 0.03961 

- 0.02054 

- 0.00199 

0.01598 

0.03346 

0.05037 

0.06683 

0.08275 

0 . 0982 : 

0.11322 

0.12781 

0.14191 

0.15547 

0.16812 




- 0.15033 

- 0.12553 

- 0.10575 

- 0.08594 

- 0.06705 

- 0.04864 

- 0.03057 

- 0 . 0128 ! 

0.00442 

0.02140 

0.03806 

0.05435 

0.07034 

0.08590 

0.10119 

0 . 1161 : 

0.13073 

0.14494 

0.15867 

0.17184 

0.18402 




- 0.09282 

- 0.07252 

- 0.05588 

- 0.03926 

- 0.02329 

- 0.00766 

0.00771 

0 . 0227 ! 

0.03766 

0.05229 

0.06671 

0.08089 

0.09482 

0.10855 

0.12203 

0.13521 

0.14816 

0.16076 

0.17300 

0.18478 

0.19570 




- 0.22623 

- 0.25766 

- 0.25881 

- 0.25715 

- 0.25046 

- 0.24091 

- 0.22859 

- 0 . 2118 ! 

- 0.19041 

- 0.16474 

- 0.13510 

- 0.10277 

- 0.07418 

- 0.04316 

- 0.01028 

0 . 0242 ! 

0.05918 

0.09269 

0.12376 

0.15196 

0.17801 




- 0.22180 

- 0.23798 

- 0.22938 

- 0.22032 

- 0.20653 

- 0.19103 

- 0.17463 

- 0 . 1569 * 

- 0.13762 

- 0.11651 

- 0.09377 

- 0.06880 

- 0.04048 

- 0.01104 

0.01729 

0.04301 

0.06803 

0.09214 

0.11529 

0.13745 

0.15910 




- 0.21743 

- 0.22212 

- 0.20847 

- 0.19620 

- 0.18074 

- 0.16395 

- 0.14652 

- 0.12811 

- 0.10906 

- 0.08898 

- 0.06796 

- 0.04625 

- 0.02392 

- 0.00143 

0.02083 

0.04271 

0.06427 

0.08517 

0.10544 

0.12510 

0.14430 




- 0.22009 

- 0.21557 

- 0.19920 

- 0.18427 

- 0.16738 

- 0.14978 

- 0.13197 

- 0.11371 

- 0.09505 

- 0.07598 

- 0.05637 

- 0.03634 

- 0.01600 

0.00444 

0.02471 

0.04481 

0.06456 

0.08392 

0.10277 

0.12103 

0.13890 




- 0.22349 

- 0.21227 

- 0.19376 

- 0.17620 

- 0.15761 

- 0.13876 

- 0.11988 

- 0.10081 

- 0.08176 

- 0.06255 

- 0.04321 

- 0.02382 

- 0.00435 

0.01505 

0.03434 

0 . 0535 ! 

0.07264 

0.09166 

0.11059 

0.12936 

0.14805 




- 0.19509 

- 0.18096 

- 0.16248 

- 0.14470 

- 0.12627 

- 0.10775 

- 0.08929 

- 0 . 0707 ' 

- 0.05226 

- 0.03375 

- 0.01518 

0.00331 

0.02187 

0.04033 

0.05875 

0 . 0771 ! 

0.09554 

0.11384 

0.13209 

0.15024 

0.16813 




- 0.15032 

- 0.13572 

- 0.11844 

- 0.10164 

- 0.08457 

- 0.06750 

- 0.05050 

- 0 . 0335 ! 


Figure A.l: Sample input data file for grid generation (continued) 


- 0.01658 

0.11762 

- 0.09282 

0.02374 

0.13894 
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0.00034 0.01715 0.03399 
0.13434 0.15109 0.16771 
0.07935 - 0.06428 - 0.04962 
0.03826 0.05271 0.06713 
0.15326 0.16758 0.18182 


0.05071 0.06742 0.08415 

0.18402 

- 0.03484 - 0.02008 - 0.00543 

0.08149 0.09586 0.11021 

0.19570 


0.10087 

0.00917 

0.12457 


Figure A.l: Sample input data file for grid generation (continued) 
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ITHETA 

IWRITE 

NCOWL 


JMAX 

IORTH 


IDEBUG 


NINPLS 

NEXPLS 

NBLPTZ 


Theta input format for blade definitions: 

if = 1, a mean camber line and tangential thickness are to be input, 
if = 2, surface theta coordinates are to be input. 

Mesh write to disk format (not used in this version). 

Cowl geometry parameter: 

if = 0 no cowl, no cowl geometry is read, 

if = 1 generate an H-grid for the ducted geometry (cowl geometry is 
read), 

if = 2 generate a multiple-block C-grid for the ducted geometry (cowl 
geometry is read). 

Defines the number of grid lines from the blade to the outer boundary 
of the axisymmetric projection of the C-grid block. 

(must be < PARAMETER KX). 

C-grid orthogonality control parameter: 

if = 0, cowl surface points for a C-grid are based on arc length inter- 
polation, 

if = 1, cowl surface points for a C-grid are adjusted for orthogonality. 

Debugging parameter: 

if = 0 no extra debugging plots, 

if = 1 several extra PLOT3D data sets are produced for debugging. 
Number of additional constant spacing inlet planes added to grid. 
Number of additional constant spacing exit planes added to grid. 
Number of points on the airfoil in the axial direction 
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NBLPTR 

NBLPTT 

NCITER 

ISOLVE 

ZINLET 

ZEXIT 

CFACLE 

CFACTE 

DFACLE 


(this must be an odd number). 

Number of points on the airfoil in the radial direction 
(this must be an odd number). 

Number of points between airfoils in the circumferential 
direction (this must be an odd number < PARAMETER KX). 
Number of iterations for the C-grid solver, 
if = 0 C-grid is generated using the algebraic solver. 

C-grid generation algorithm control parameter: 
if = 0 C-grid is generated using the Poisson solver, 
if = 1 C-grid is generated using the variational solver. 

Axial location of the initial inlet boundary plane 
(before constant spacing inlet planes are added). 

Axial location of the initial exit boundary plane 
(before constant spacing exit planes are added). 

Factor which determines the axial extent of the C-grid 
upstream of the cowl leading edge (multiplied by 
cowl axial chord). 

Factor which determines the axial extent of the C-grid 
downstream of the cowl trailing edge (multiplied by 
cowl axial chord). 

Factor which determines the radial width of the C-grid 
upstream of the cowl leading edge (if > 1, then 
larger than tip gap, if < 1, then smaller than tip 
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DFACTE 

RATIN 

RATEX 

RATBB 

RATBLZ 

RATBLR 

RATBLT 

RATTOB 

NHUB 

NTIP 


gap). 

Factor which determines the radial width of the C-grid 
downstream of the cowl trailing edge (if > 1, then 
larger than tip gap, if < 1, then smaller than tip 

gap). 

Maximum axial ratio of cells in the inlet region 
(must be > 1). 

Maximum axial ratio of cells in the exit region 
(must be > 1). 

Maximum axial ratio of cells between blade rows 
(must be >1). 

Maximum axial ratio of cells on blades 
(must be > 1). 

Maximum radial ratio of cells on blades 
(must be > 1). 

Maximum circumferential ratio of cells blade to blade 
(must be > 1). 

Maximum radial ratio of cells from tip (unducted H-grids) or 
the cowl upper surface (ducted H-grids) or the C-grid upper 
boundary (ducted C-grids) to the outer boundary 
(must be > 1). 

Number of X,R hub coordinate input pairs to be read. 
Number of X,R outer boundary coordinate pairs to be read. 



113 


XHUB,RHUB Hub coordinates. 

All of the hub axial coordinates are given first (8 per 
line), followed separately by the radial coordinates. 
The number of points input is determined by NHUB. 
XTIP,RTIP Outer boundary coordinates. 

All of the outer axial coordinates are given first (8 per 
line, followed separately by the radial coordinates. 

The number of points input is determined by NTIP. 


The following 5 variables are used only for a ducted geometry. They should be 
deleted for an unducted case. 


NSHR 

NSHRIN 


NSHROU 


NPBCAB 

XSHR,RSHR 


Number of X,R cowl coordinate input pairs to be read. 
Number of axial grid lines between first blade row 
leading edge and cowl leading edge. 

(this must be an odd number). 

Number of axial grid lines between last blade row 
trailing edge and cowl trailing edge. 

(this must be an odd number). 

Number of radial grid lines between lower cowl surface 
and blade tip. 

Cowl coordinates. The coordinates are input in a 
wraparound fashion beginning at the trailing edge 
and proceeding in the clockwise direction. The 
trailing edge point must be duplicated as the last coordinate. 
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All of the cowl axial coordinates are given first (8 per 
line), followed separately by the radial coordinates. 
The number of points input is determined by NSHR. 


The remaining variables define the propfan blades themselves. The blades are 
defined by radial slices with identical axial and radial coordinates on each side of the 
airfoil. The slices are given from the hub to the tip. The tip cross section must be 
specified in a manner which provides some clearance between the tip of the blade 
and the cowl lower surface, even if the two parts are actually mated. The 8 variables 
which follow are duplicated for each blade row when more than one row is considered. 


NBLA 

NBLCRC 

NPPRC 

IPCH 


ZPCA 


THPCA 


Number of blades in this blade row. 

Number of radial cross section used to define the blade. 
Number of axial points per radial cross section used to 
define the blade. 

Pitch change trigger: 

if = 0 no pitch change is performed, 

if = 1 the pitch change variables ZPCA, THPCA, and BETA 
are read in and the blade setting angle is adjusted 
to the value specified by BETA. 

Axial location of the blade pitch change axis. 

The pitch change axis is assumed perpendicular to 
the axis of rotation. 

Circumferential position of the blade pitch change axis. 

The pitch change axis is assumed perpendicular to 
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the axis of rotation. 

BETA Absolute 3/4 radius blade setting angle. (Measured positive 
from the plane of rotation. The blade setting angle reference 
plane is the same for either positive or negative rotation 
. [see Fig. A.2]). 

ZBLA,RBLA, 

THBLA,TTBLA Blade coordinates. 

Each of the blade axial coordinates is given first 
(8 per line) for each cross section, cross sections are 
given from hub to tip) followed by all the radial 
coordinates. The blade circumferential coordinates are given 
next, and may be in two different formats depending 
on the variable ITHETA. The preferred format is 
ITHETA=2, where the actual tangential surface 
coordinates are given for the surfaces ordered 
in the counter-clockwise direction (for a 
positive rotation, this would be the suction 
surface first, then the pressure surface). 

The grid generation system produces two primary types of output files: the 
standard output and the grid or grids (depending on whether a C-grid or H-grid is 
being produced) themselves. A description of these files is given in Table A.l. 

A typical H-grid generated from the sample data set given in Fig. A.l is shown 
in Fig. A. 3. The blade tip gap has been expanded to illustrate the grid in this region. 
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Figure A. 2: Blade setting angle reference plane description 
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Table A.l: Description of input/output files and UNIX-based filenames for H-grid 
generation 


Name 

Lu No. 

Description 

fort .4 

4 

Unformatted grid data set for first blade row. 
This may be viewed using PL0T3D, or used 
directly as input for the aerodynamic 
analysis. 

fort. 14 

14 

Unformatted grid data set for second blade row. 

fort .24 

24 

Unformatted grid data set for third blade row (etc.). 

fort .98 

98 

Geometry and input for grid generator. 

fort.99 

99 

Standard output from grid generator. 


The C-grid generation scheme utilizes the same input and standard output files, 
but produces a number of additional files which are described in Table A.2. 

A typical C-grid generated from the sample input data file given in Fig. A.l 
modified by the condition NCOWL = 2 is shown in Fig. A.4. 


A. 3 Subroutine Description 

A list of the grid generation program subroutines and their functions is given 
below for reference. A skeleton program flowchart is illustrated in Fig. A. 5. 

SUBROUTINE DESCRIPTION 


CHGRID 

ALPHM 

CCOWL 

CGRID 


Main calling routine. 

Plotting routine - specifically for alphanumeric labels, etc. 
Routine for setting up and constructing C-grid boundaries. 
Routine for determining C-grid interior point distribution. 


Figure A. 3: Sample grid based on H-grid generation 
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Figure A.4: Sample grid based on multiple- block C-grid generation 
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Table A.2: Description of input/output files and UNIX-based filenames for C-grid 



generation 


Name 

Lu No. 

Description 

fort. 4 

4 

Unformatted grid data set of axisymmetric projection 
of coupled grid blocks without C-grid. 

fort. 15 

41 

PostScript plot data of C-grid quality. 

fort .41 

41 

Block #1 unformatted grid data set. 

fort .42 

42 

Block #2 unformatted grid data set. 

fort .43 

43 

Block #3 unformatted grid data set. 

fort .44 

44 

Block #4 unformatted grid data set. 

fort .45 

45 

Block #5 unformatted grid data set. 

fort .98 

98 

Geometry and input for grid generator. 

fort. 99 

99 

Standard output from grid generator. 

fort. 17 

17 

Unformatted grid data set for debugging. 

fort.27 

27 

Unformatted grid data set for debugging. 

fort. 37 

37 

Unformatted grid data set for debugging. 

fort. 47 

47 

Unformatted grid data set for debugging. 

fort. 57 

57 

Unformatted grid data set for debugging. 

fort .67 

67 

Unformatted grid data set for debugging. 
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COWEX 

COWIN 

CURPLT 

END 

ERROR 

EULER 

GENBB 

GENBL 

GENIN 

GUESS 

INRSCT 

INTERP 

METRIC 

NUMBER 

NUMPTS 

POISON 

POLFIT 

PRINTO 

SCALE 


Routine which determines the outer boundary of the exit region. 
Routine which determines the outer boundary of the inlet region. 
Plotting routine - specifically for plotting data as curves. 

Plotting routine - specifically for ending a plot. 

Routine for internal evaluation of potential input data 
or grid errors. 

Variational-based interior point solver for C-grids. 

Routine for generating grid point distribution between blades. 
Routine for determining blade grid point distribution. 

Routine for determining inlet region grid point distribution. 
Routine for generating an initial guess for C-grid interior 
point grid generation scheme. 

Spline-based routine for determining the intersection of two 
curves. 

Polynomial-based interpolation routine. 

Routine for determining metric terms in C-grid generation. 
Plotting routine - specifically for plotting numbers. 

Routine for determining the number of grid lines required for 
a specified grid cell ratio in a given region. 

Poisson-based interior point solver for C-grid. 

Polynomial-based interpolation routine. 

Output routine for grid generation. 

Plotting routine - specifically for generating scales. 
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SPLINT Cubic spline polynomial interpolation routine. 

START Plotting routine - specifically for initiating a plot. 

WEIGHT Grid point weighting routine for variational grid generation. 
ZCOWL Routine for generating H-grids about cowls. 



CHGRID 

I 

— PRINT 123 

I 

— ERROR 
I 

— INTERP 
I 

— GENBL 

— COWIN 
I 

— GENIN 
I 

— GENBB 
I 

— COWEX 
I 

— GENEX 
I 

— ZCOWL 
I I 

I — COWIN 

I I 

I — COWEX 

I I 

I — CCOWL 

I I 

I — CGRID 

I I 

I — GUESS 

I I 

I — WEIGHT 

I I 

I — METRIC 

I I 

I — POISON 

I I 

I — EULER 

I I 

I — ANALYZ 

I I 

I — START 

I I 

I — SCALE 

I I 

I — CURPLT 

I — ALPHM 

I I 

I — END 


— NUMPTS 
I 

— GENBL 
I 

— GENIN 
I 

— GENBB 
I 

— GENEX 
I 

STOP 


Figure A.5: Program flowchart for ducted propfan analysis grid generation 
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APPENDIX B. 3D EULER SOLVER OPERATING INSTRUCTIONS 


B.l Introduction 

This appendix contains the computer program User’s Manual for the 3D Euler 
ducted propfan aerodynamic analysis. Again, as in the grid generation, both an H- 
grid and a multiple-block C-grid solution scheme have been developed for the analysis 
of ducted propfans. In this case, separate programs are used for each of these analyses, 
although a number of similarities exist for the operation of each code. The H-grid 
flow solver can also be used as is for the analysis of unducted propfans. The H-grid 
version also has the capability of analyzing counter-rotating ducted propfans using 
the average-passage system of equations. 

The flow solver source programs are written in FORTRAN 77, and have been 
used successfully on Cray UNICOS and IBM VM/CMS mainframe computer systems 
as well as Silicon Graphics 4D workstations using a UNIX operating system. Array 
dimensions in the programs are specified by a PARAMETER statement in each code. 
The H-grid analysis parameters are specified in the statement: 

PARAMETER( IMX=125, JMX=25, KMX=16 ) 
which appears in every subroutine. 

The H-grid Euler analysis PARAMETER variables are defined as: 
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IMX Number of axial grid lines + 1 
JMX Number of radial grid lines + 1 
KMX Number of circumferential grid lines + 1 

The C-grid Euler analysis code uses several PARAMETER statements all in the 
main routine given as: 

PARAMETER( NBLKS = 5 ) 

PARAMETER( IMX1=125, JMX1=25, KMX1=16 ) 
PARAMETER( IMX2=125, JMX2=25, KMX2=16 ) 
PARAMETER( IMX3=125, JMX3=25, KMX3=16 ) 
PARAMETER( IMX4=125, JMX4=25, KMX4=16 ) 
PARAMETER( IMX5=125, JMX5=25, KMX5=16 ) 

The C-grid Euler analysis PARAMETER variables are defined as: 

NBLKS Number of grid blocks (this need not be changed) 

IMX1 Number of axial grid lines in block #1 + 1 
IMX2 Number of axial grid lines in block #2 + 1 
IMX3 Number of axial grid lines in block #3 + 1 
IMX4 Number of axial grid lines in block #4 + 1 
IMX5 Number of axial grid lines in block #5 + 1 
JMX1 Number of radial grid lines in block #1 + 1 
JMX2 Number of radial grid lines in block #2+1 
JMX3 Number of radial grid lines in block #3+1 
JMX4 Number of radial grid lines in block #4+1 
JMX5 Number of radial grid lines in block #5+1 
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KMX1 Number of circumferential grid lines in block $1 + 1 

KMX2 Number of circumferential grid lines in block $2 + 1 

KMX3 Number of circumferential grid lines in block $3 + 1 

KMX4 Number of circumferential grid lines in block $4 + 1 

KMX5 Number of circumferential grid lines in block $5 + 1 

The requirement that the PARAMETER variables be 1 element larger than the 
grid dimensions results from the use of phantom points outside the computational 
domain to impose the numerical boundary conditions. 

An important requirement here is that the first array dimension must be larger 
than the second and third array dimensions. This is based on the dimensioning scheme 
used for certain temporary variables in the code. This condition therefore imposes 
some limitations on the grid sizes that may be run without modifying the code. Under 
certain conditions, it may also be required that the second grid dimension be larger 
than the third. 

Approximate computational storage and CPU requirements can be estimated for 
either the E-grid or C-grid Euler codes from the following formulas: 

CPU sec ss 4.0xl0~ **($ grid points)($ iterations) 

Memory Mw « 2.7xl0“'*($ grid points) 

These formulas are valid for a Cray-II computer operating under the UNICOS 
environment and the cf77 compiler. 

The programs produce plotting output suitable for viewing using the PLOT3D 
graphics software package developed at the NASA Ames Research Center. PLOT3D 
output is written to FORTRAN unit 9 for the H-grid solver, while the C-grid version 
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writes individual block grid information to units 91*95 for the five-block architecture. 

B.2 Input/Output Files 

An example input data file for the H-grid and multiple-block C-grid 3D Euler 
solvers is shown in Fig. B.l. A brief description of the form of the data file follows. 
It should be noted that the H-grid Euler solver was developed from an existing 3D 
Euler code with a number of acceleration features including axisymmetric averaging, 
reduced solution domain, grid coarsening and grid refinement. Many of the input 
variables listed below deal with controlling these functions. Due to the nature of the 
analysis for ducted propfans, however, several of these algorithms are not enabled. 
Therefore, a description of the controlling variables is provided for reference only, and 
should not be considered functional in the present analysis. 

VARIABLE DESCRIPTION 

MACH Freestream Mach number. 

GAMMA Specific heat ratio (cp/cv). 

PEXIT Freestream static pressure. This is applied at the 
outer exit boundary and integrated radially inward 
for propellers (DUCT=0.0), or applied at the hub and 
integrated radially outward for compressors (DUCT=1.0) 
to satisfy radial equilibrium. 

OMEGA Nondimensional rotational speed (this value will be 

adjusted to match the advance ratio when ADVR ^ 0.0). 
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+ — MACH — ++-GAMMA — ++-PEXIT — ++-OMEGA — ++ — ADVR — ++ — DUCT — ++-COWL + 

0.8000 1.4000 0.65602 0.0000 -3.0600 0.0 1.0 

+--CFL ++ — VIS2 — ++ — VIS4 — ++ — QFIL — ++ BC ++-BLDROW-++-DTAXI — + 

-5.5 2.00 1.000 0.0000 0.0000 1.0 0.80 

+ -FNCMAX-++ — REST— ++— SAVE — ++-FISTEP-++-FNPRNT-++-ROWMAX — H+-FPRAXI-+ 

720.0 1.0 1.0 1.0 1.0 1.0 0.0 

H — REFINI-+H — REFINJ-+H — REFINK — — EPSX — ++ — EPSY — ++ — EPSZ — ++-FPRBF— + 

1.0 1.0 1.0 2.00 2.00 2.00 0.0 

+-RFINII-++-RFINJJ-++-RFINKK-++ — FAXI — ++-FAXIMX-++-FR3DMX-++ — FTURB-+ 

1.0 1.0 1.0 0.0 99999.0 0.0 0020.0 

+-FAXIS — ++-FAXIE — ++-FSTART-++-FEXIT — ++-JBASE — ++-OMEGS1-++-OMEGS2-+ 
-10.0 -9.0 5.00 4.00 1.0 0.0000 0.0000 


Figure B.l: Sample input data file for ducted propfan 3D Euler solver 
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ADVR Advance ratio (J = U/nD): 

if = 0.0, rotational 6peed is determined by OMEGA. 

DUCT Internal flow duct option: 

if = 0.0, external flow options are utilized (propeller), 
if = 1.0, an internal flow is assumed (compressor). No slip 
boundary conditions are applied at the outer 
boundary rather than the characteristic far field 
condition. 

COWL Cowl geometry trigger: 

if = 0.0, an unducted geometry is assumed, 
if = 1.0, duct boundary conditions are applied 
(only meaningful for H-grid flow solver). 

CFL Time step parameter: 

if < 0 then this is the CFL number used 
to determine the calculation time step for 
local time stepping, 

if > 0.0 then this is the CFL number used 
to determine the calculation time step without local 
time stepping. 

VIS2 Second order damping coefficient (~ 2.0) 

(divided by 4 in the code) 

VIS4 Fourth order damping coefficient (w 1.0) 

(divided by 64 in the code) 
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QFIL 

BC 


BLDROW 


DTAXI 

FNCMAX 

REST 


SAVE 


FISTEP 

FNPRNT 

ROWMAX 


Not active in this version (0.0). 

Exit boundary condition trigger: 

if = 0.0, characteristic boundary condition based on is used, 
if = 1.0, characteristic boundary condition based on mass flow is used. 
Blade row parameter. This value determines which blade row 
is being calculated during a multiple blade row solution. 

For single blade row calculations, this should be = 1.0. 

Not active in this version. 

Maximum number of time steps to be performed. 

Restart option parameter: 

if = —1.0 a restart file and body force file are read, 
if = 0.0 no restart file, initialize variables in code, 
if = 1.0 a restart file is read. 

Save restart file option parameter: 

if = —1.0 restart files and body source terms are output 
at the end of a run, 

if = 0.0 no restart file is output at the end of the run, 
if = 1.0 a restart file is output at the end of the run. 

Not active in this version (= 1.0). 

Not active in this version (= 0.0). 

Maximum number of blade rows in the current solution. 

For single rotation, = 1.0. For multiple blade rows, 
this parameter indicates that additional grids and 
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FPRAXI 


REFINI 

REFINJ 

REFINK 

EPSX 

EPSY 

EPSZ 

FPRBF 


RFINII 

RFINJJ 

RFINKK 

FAXI 

FAXIMX 

FR3DMX 


body source terms are to be read. 

Print axisymmetric solution, 

if = 0.0 axisymmetric solution is not printed, 

if = 1.0 axisymmetric solution is printed. 

Grid-coarsening trigger, not active in this version (1.0). 
Grid-coarsening trigger, not active in this version (1.0). 
Grid-coarsening trigger, not active in this version (1.0). 

Implicit residual smoothing coefficient in the axial 
direction (« 2.0 is a typical value). 

Implicit residual smoothing coefficient in the radial 
direction (ss 2.0 is a typical value). 

Implicit residual smoothing coefficient in the circumferential 
direction (w 2.0 is a typical value). 

Print trigger for body forces: 

if = 0.0, body forces are not printed, 

if = 1.0, body forces are printed. 

Fine mesh conversion, not active in this version (1.0). 

Fine mesh conversion, not active in this version (1.0). 

Fine mesh conversion, not active in this version (1.0). 

Number of axisymmetric solutions, not active in this version (0.0). 
Number of 3D passes before axisymmetric solutions, not active in this 
version (99999.0). 

Number of reduced domain 3D passes, not active in this version (0.0). 
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FTURB Not active in this version (0.0). 

FAXIS Not active in this version (0.0). 

FAXIE Not active in this version (0.0). 

FSTART Not active in this version (0.0). 

FEXIT Not active in this version (0.0). 

JBASE Not active in this version (0.0). 

OMEGS1 Not active in this version (0.0). 

OMEGS2 Not active in this version (0.0). 

The Euler solver system produces three primary types of output files: the stan- 
dard output, plot output, and restart files. A list and description of each of these 
files is given in Table B.l for the H-grid Euler analysis and Table B.2 for the C-grid 
Euler analysis. 

B.3 H-grid Euler Solver Subroutine Description 

A list of the 3D H-grid Euler solver program subroutines and their functions is 
given below for reference. A skeleton program flowchart is illustrated in Fig. B.2. 

SUBROUTINE DESCRIPTION 

HPR03D Main calling routine. 

BAXI Calculates axisymmetric flow quantities. 

BAXI2D Calculates axisymmetric flow quantities for 2D solution. 

BAXIA Axisymmetric flow averaging routine. 

BCCL Boundary condition routine for cowl lower surface. 
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Table B.l: Description of input/output files and UNIX-based filenames for H-grid 

Euler solver 



Lu No. 

Description 

fort.l 

1 

New restart file output for first blade row. 

fort .2 

2 

Old restart file input for first blade row. 

fort. 3 

3 

Old body force output for first blade row. 

fort .4 

4 

Unformatted grid data set for first blade row. 

This may be viewed using PL0T3D, or used directly 
as input for the aerodynamic analysis. 

fort. 5 

5 

Main program input data file (may be redirected). 

fort. 6 

6 

Main program output data file (may be redirected). 

fort. 9 

9 

Unformatted flow data set. This may be viewed 
using PL0T3D. 

fort. 11 

11 

New body force output for current blade row. 

fort. 12 

12 

Old restart file for second blade row. 

fort.22 

22 

Old restart file for third blade row (etc.). 

fort. 13 

13 

Old body force output for second blade row. 

fort.23 

23 

Old body force output for third blade row (etc.). 

fort. 14 

14 

Unformatted grid data set for second blade row. 

fort. 24 

24 

Unformatted grid data set for third blade row (etc.). 
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Table B.2: Description of input /output files and UNIX-based filenames for multiple- 
block C-grid Euler solver 


Name 

Lu No. 

Description 

fort. 11 

11 

New restart file output for block $1. 

fort. 12 

12 

New restart file output for block #2. 

fort. 13 

13 

New restart file output for block #3. 

fort. 14 

14 

New restart file output for block #4. 

fort. 15 

15 

New restart file output for block #5. 

fort. 21 

21 

Old restart file input for block #1. 

fort. 22 

22 

Old restart file input for block #2. 

fort. 23 

23 

Old restart file input for block #3. 

fort. 24 

24 

Old restart file input for block #4. 

fort. 25 

25 

Old restart file input for block #5. 

fort.41 

41 

Block #1 unformatted grid data set. 

fort. 42 

42 

Block #2 unformatted grid data set. 

fort .43 

43 

Block #3 unformatted grid data set. 

fort. 44 

44 

Block #4 unformatted grid data set. 

fort .45 

45 

Block #5 unformatted grid data set. 

fort. 51 

51 

Main program input data file for block #1. 

fort. 52 

52 

Main program input data file for block #2. 

fort. 53 

53 

Main program input data file for block #3. 

fort. 54 

54 

Main program input data file for block #4. 

fort. 55 

55 

Main program input data file for block #5. 

fort .6 

6 

Main program output data file (may be redirected). 

fort. 61 

61 

Block #1 printed output file. 

fort. 62 

62 

Block #2 printed output file. 

fort .63 

63 

Block #3 printed output file. 

fort. 64 

64 

Block #4 printed output file. 

fort .65 

65 

Block #5 printed output file. 

fort.91 

91 

Unformatted flow plot data set for block #1. 

fort. 92 

92 

Unformatted flow plot data set for block #2. 

fort. 93 

93 

Unformatted flow plot data set for block #3. 

fort. 94 

94 

Unformatted flow plot data set for block #4. 

fort. 95 

95 

Unformatted flow plot data set for block #5. 

All unformatted grid and flow plot data sets may be 
viewed using PL0T3. 
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BCCU 

BCEXT 

BCINL 

BCJ 

BCJ2 

BCK 

BCWAL2 

BCWALL 

BFORCE 

BFREAD 

ERROR 

ERROR2 

ERROR3 

ERRORP 

EULER 


Boundary condition routine for cowl upper surface. 

Boundary condition routine for exit cells. 

Boundary condition routine for inlet cells. 

Boundary condition routine for j=constant surfaces (hub and 
outer boundaries). 

Boundary condition routine for j=constant surfaces (hub and 
outer boundaries for 2D solution). 

Boundary condition routine for the blade surfaces. 
Axisymmetric boundary conditions at the exit for the 
2D solution. 

Boundary condition routine for periodic cells. 

Routine for calculating body force terms using the 
axisymmetric equations. 

Routine to read in body force terms for all adjacent blade 
rows. 

Convergence-checking routine for Runge-Kutta solver (not used 
for propfans). 

Convergence-checking routine for 2D Runge-Kutta solver. 
Convergence-checking routine for reduced domain 
3D solution. 

Convergence-checking routine for Runge-Kutta solver for 
propfan calculations. 

Runge-Kutta solver. 
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EULER2 

EULER3 

FILTE2 

FILTE3 

FILTER 

FIXUP 

GRIDG 

INITA 

INITB 

INPUT 

METRI2 

METRIC 

0P2D 

OUTPUT 

PTCALC 

REFIN 

REFIN2 

REFIN3 

RESID 

RESID2 


2D Runge-Kutta solver. 

Reduced domain Runge-Kutta solver. 

2D artificial dissipation routine. 

Reduced domain artificial dissipation routine. 

Artificial dissipation routine. 

Routine to correct boundary cells before output averaging 
is performed. 

Routine to read and set up the grid. 

Routine to initialize flowfield with axisymmetric solution. 
Routine to replace flowfield with latest axisymmetric solution. 
Routine to read in input data and set up reference values. 
Routine to calculate cell volumes and surface normals for 
2D solution. 

Routine to calculate cell volumes and surface normals. 

Routine to calculate axisymmetric residual of axisymmetric 
flow variables. 

Routine to print output and save restart files. 

Routine to determine inlet total pressure profile. 

Routine to refine a 3D solution. 

Routine to refine a 2D solution. 

Routine to refine a reduced domain 3D solution. 

Implicit residual smoothing routine. 

2D implicit residual smoothing routine. 
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RUNGE 

RUNGE2 

STEP 

STEP2 

TEST2D 

VOLUX 

VOLUX2 


Runge-Kutta flux calculation routine. 

2D Runge-Kutta flux calculation routine. 

Routine to determine time step. 

Routine to determine time step for 2D solution. 
Routine to calculate the axisymmetric average of the 
residuals of the 3D flowfleld. 

Routine to calculate an individual cell volume. 
Routine to calculate an individual 2D cell volume. 


B.4 Multiple-Block C-Grid Euler Solver Subroutine Description 

A list of the multiple-block C-grid 3D Euler solver program subroutines and their 
functions is given below for reference. A skeleton program flowchart is illustrated in 
Fig. B.3 for the multiple-block C-grid 3D Euler solvers. 

SUBROUTINE DESCRIPTION 


CPR03D 

BCI12 

BCI13 

BCI14 

BCI23 

BCI25 

BCI33 

BCI34 


Main calling routine. 

Boundary condition routine between block 1 and 2. 
Boundary condition routine between blocks 1 and 3. 
Boundary condition routine between blocks 1 and 4. 
Boundary condition routine between blocks 2 and 3. 
Boundary condition routine between blocks 2 and 5. 
Boundary condition routine along C-grid branch cut. 
Boundary condition routine between blocks 3 and 4. 



HPR03D 


— INPUT 
I I 

I — GRIDG 

I I 

I — PTCALC 

I I 

I — INIT 

I I 

I — INITA 

I I 

I — BFREAD 

I I 

I — INITB 

I I 

I — METRIC 

I I I 

I I — VOLUX 

I I 

I — METRI2 

I I 

I — VOLU2 

— STEP 
I 

-- EULER 
I I 

I — FILTER 

I I 

I — RUNGE 

I I I 

I I — RESID 

I I 

I — BCINL 

I I 

I — BCWALL 

I I 

I — BCEXT 

I I 

I — BCJ 

I I 

I — BCK 

I I 

I — BCCL 

I I 

I — BCCU 

I I 

I — FIXUP 

I 

— ERRORP 
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— EULER2 

— BAXI2D 
I 

— TEST2D 
I 

— OP2D 
I I 

I — BCJ2 
I I 

I — BCWAL2 
I 

— STEP2 
I 

— FILTE2 
I 

— RUNGE2 
I I 

I — RESID 

I 

— BCWAL2 
I 

— BCJ2 

— ERROR2 
I 

— BCINL 
I 

— BCEXT 
I 

— BCJ 
I 

— BCK 


EULER3 
I 

— STEP 
I 

-- FILTER 

— RUNGE 
I I 

I — RESID 

I 

— BCINL 
I 

— BCWALL 
I 

— BCEXT 
I 

— BCJ 
I 

— BCK 
I 

— ERROR3 

— REFIN 
I I 

I — GRIDG 

I I 

I -- REFIN3 

I I 

I — REFIN2 

I I 

I — METRIC 

I I 

I — VOLUX 

I 

— OUTPUT 
I 

— BAXI 
I 

— TEST2D 
I 

— BFORCE 
I I 

I -- BCJ2 

I I 

I — BCWAL2 

I 

STOP 


I 


Figure B.2: Program flowchart for ducted propfan analysis H-grid Euler solver 
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BCI35 

BCI45 

BCINL 

BCJ 

BCK 

BCPER 

ERRORP 

EULER 

FILTER 

FIXUP 

GRIDG 

INIT 

INPUT 

METRIC 

OUTPUT 

PTCALC 

REFIN 

RESID 

RUNGE 

STEP 


Boundary condition routine between blocks 3 and 5. 
Boundary condition routine between blocks 4 and 5. 
Boundary condition routine for inlet cells. 

Boundary condition routine for j=constant surfaces (hub and 
outer boundaries). 

Boundary condition routine for the blade surfaces. 

Boundary condition routine for periodic cells. 
Convergence-checking routine for Runge-Kutta solver for 
propfan calculations. 

Runge-Kutta solver. 

Artificial dissipation routine. 

Routine to correct boundary cells before output averaging 
is performed. 

Routine to read and set up the grid. 

Routine to initialize solution. 

Routine to read in input data and set up reference values. 
Routine to calculate cell volumes and surface normals. 
Routine to print output and save restart files. 

Routine to determine inlet total pressure profile. 

Routine to refine a 3D solution. 

Implicit residual smoothing routine. 

Runge-Kutta flux calculation routine. 

Routine to determine time step. 
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VOLUX Routine to calculate an individual cell volume. 
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CPR03D 

I 


— INPUT 
1 1 

i 

i 

i 

— BCI33 

| — GRIDG 

1 1 

i 

i 

1 

— ERRORP 

| — PTCALC 

1 1 

i 

i 

1 

— FIXUP 

1 — INIT 

1 1 

i 

i 

1 

— OUTPUT 

1 — METRIC 

1 1 

i 

i 

1 

— P3DOUT 

1 VOLUX 

1 

i 

i 

1 

STOP 

— STEP 

i 


I 

-- FILTER 

i 

1 

1 

i 


1 

-- RUNGE 

i i 

1 

1 

i 


1 1 

| — RESID 

i 

1 

1 


i 

— BCPER 

— BCBLA ' 

i 

i 

i 

i 


1 

— BCCOWL 
1 

1 

1 

i 


1 

-- BCHUB 
1 

1 

1 


1 

— BCINL 
| 

1 

1 


1 

— BCEXT 

i 

1 

i 

i 


1 

— BCOUT 

i 

1 

i 


1 

-- BCI12 

i 

i 

i 

i 


I 

-- BCI13 

1 

i 

i 


-- BCI14 

1 

i 

i 


1 

~ BCI23 

I 

1 

i 

i 


— BCI34 

i 

1 

i 

i 


1 

-- BCI45 

I 

i 

i 


1 

-- BCI25 
1 

1 

i 

i 



Figure B.3: Program flowchart for ducted propfan analysis multiple-block C-grid 

Elder solver 
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